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Abstract. This article investigates the full Boltzmann equation up to second 
order in the cosmological perturbations. Describing the distribution of polarized 
radiation by a tensor valued distribution function, we study the gauge dependence 
of the distribution function and summarize the construction of the gauge-invariant 
distribution function. The Liouville operator which describes the free streaming 
of electrons, and the collision term which describes the scattering of photons 
on free electrons are computed up to second order. Finally, the remaining 
dependence in the direction of the photon momentum is handled by expanding 
in projected symmetric trace-free multipoles and also in the more commonly 
used normal modes components. The results obtained remain to be used for 
computing numerically the contribution in the cosmic microwave background 
bispcctrum which arises from the evolution of second order perturbations, in 
order to disentangle the primordial non-Gaussianity from the one generated by 
the subsequent non-linear evolution. 



Introduction 

The cosmic microwave background (CMB) has become in the past twenty years a 
central observable of modern cosmology. The properties of the CMB temperature 
fluctuations depend both on the initial conditions set at the end of the primordial in- 
flationary era, and on their evolution through time in the post-inflationary eras. The 
theory of cosmological perturbations around a space-time with maximally symmetric 
spatial sections is a cornerstone of our understanding of the large scale structure of the 
universe. The relativistic matter (photons, neutrinos) is described with a statistical 
approach [UEIE] , also referred to as kinetic theory, in which we use a distribution func- 
tion whose evolution is given by a Boltzmann equation. As for non-relativistic matter 
[baryons, cold dark matter], a fluid approximation is usually sufficient. Given the 
typical amplitude (10~^) of the metric fluctuations, the dynamics of the cosmological 
perturbations was so far mainly studied at linear order around a Friedmann-Lemaitre 
(FL) space-time, and characterized statistically by the power spectrum. This method 
allows to relate the CMB angular power spectrum to the initial power spectrum, which 
opens a window on the early universe. However, the linear order of perturbations fails 
to capture the intrinsic non-linear features of General Relativity which enter both 
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the initial conditions and the evolution. The CMB measurements, though already of 
high precision [H [5] , are soon going to improve with the new forthcoming missions 
such as Planck [S], and will be sensitive to these non-linear effects. It becomes thus 
necessary to go beyond this linear perturbation scheme by studying the second-order 
perturbations. This is even more crucial if we want to estimate the bispectrum in 
the cosmic microwave background since this can only arise from non-Gaussian initial 
conditions set by inflation or from non-linear evolution. In order to improve our the- 
ory of the early universe and discriminate between different models of inflation, it is 
thus necessary to disentangle the primordial non-Gaussianity from the one induced by 
non-linear evolution. We thus need to extend the program followed at first order in 
perturbations up to second order. 

Since cosmological perturbations are plagued by the gauge freedom, we need to 
build a full set of gauge invariant perturbation variables and derive the perturbed 
Boltzmann and Einstein equations in terms of these variables to obtain their dy- 
namical equations. The first-order gauge invariant perturbation variables in the fluid 
approach were built in [7], and in the kinetic theory in [HI [9]. In the inflationary 
era and for slow-roll one field inflation, the quantization in the linear equations of 
the canonical degrees of freedom [10] which transfer to Gaussian classical fluctuations 
enable us to fix the initial conditions for the post-inflationary dynamics on super- 
Hubble scales. By using quantities which are conserved for modes larger than the 
Hubble radius [11] [121 IE] , we can ignore the details between the end of inflation and 
the subsequent eras. In this program, the evolution of radiatiorQ requires a special 
care since the distribution function is shaped by Compton scattering by free electrons 
and in particular this generates polarization. At first order, the evolution equations 
describing polarized radiation were studied intensively in [Mj [151 [SI j and codes have 
been made available [17[ 118] for integrating numerically these equations and thus an- 
alyzing the CMB data. 

At second-order, the program remains so far incomplete and this paper aims at 
filling this gap. The quantization of the canonical degrees of freedom in the models 
of slow-roll one field inflation with non-linear couplings has been investigated in the 
interaction picture [HI [20] up to the loop corrections [211 EH aiming at computing 
the level of primordial non-Gaussianity. This generically predicts negligible amounts 
of non-Gaussianity [H] whereas multi-field inflation [24l[2ll[26l[27l[28l[29l[30l[3ll|32] 
can generate significant levels of non-Gaussianity. Since the use of adapted estima- 
tors [33l [34l [35] tends to show that there might be a detectable amount of non- 
Gaussianity in the CMB [H [5] [36] , the understanding of non- linear evolution is of first 
importance in order to constrain these models of inflation. In particular all non-linear 
effects in the foreground [371 IMl [Ml [iQ] (see [41] for a review) have to be understood 
and estimated. In order to develop the mathematical tools for this general study, the 
gauge issue in the fluid approach was studied in [42] and gauge invariant variables 
were built in [12l [43] up to second order in perturbations. This fluid approximation 
has already been used to understand the general form of the bispectrum generated by 
evolutionary effects [44]. As for the more general kinetic approach, the gauge issue was 
studied in our previous paper [J^ , and the evolution equations through free-streaming 

X Usually, the word radiation is used for rclativistic particle. In the context of CMB, it is often used 
just for photons, as it is the case in this paper. 
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and Compton collision on free electrons were derived in 133 UHl SH] , but it is so 
far restricted to unpolarized radiation, which is inconsistent since Compton scattering 
does generate polarization. In this article, we will extend these works to polarized 
radiation, leaving open the issue of second-order numerical integration. It should be 
mentioned that an alternative program consists in working directly with covariantly 
defined quantities in the so-called 1-1-3 covariant formalism. The first order canonical 
degrees of freedom required to fix the initial conditions in this approach have been 
identified [50], the conserved quantities, used to ignore the detail between the inflation- 
ary era and the subsequent eras, have been built [5T1 [5^ and the dynamical equations 
up to second order [53l [54l [55l [56l [57l [58l [59] have also been extensively studied, but 
so far leaving aside for the second order the problem of mode decomposition and the 
treatment of polarized polarization. The results presented in this paper can be used in 
a straightforward manner for the description of collisions in this 1-1-3 formalism, and 
we even correct a small mistake for the unpolarized case. The two formalisms should 
of course lead to the same conclusion and they have been compared in [501 [Ml [S2] ■ 

This paper is organized as follows. We first review the description of polarized 
radiation with a tensor valued distribution function in section [TJ We then recall 
the gauge dependence and the construction of gauge invariant variables when using 
a fluid approximation [3 In section [31 we summarize the results of our previous 
paper [45] which focused on the gauge dependence of the distribution function and the 
construction of a gauge invariant distribution function in the kinetic theory, and we 
extend the results to the tensor valued distribution function. This formalism is then 
used to compute the second order gauge invariant Liouville operator, for radiation in 
section [31 and for matter in section [S] The collision term for both cases is computed in 
section [6l Finally, in section [71 we develop the necessary tools to express these results 
in terms of the normal modes decomposition. 



1. Kinetic theory 

1.1. Momentum and tetrad 

In the kinetic description of radiation, the momentum of photons is usually 
decomposed onto an orthonormal basis, that is using a tetrad field. Though not 
compulsory, this facilitates the separation between the magnitude of the momentum 
and its direction represented by a spacelike unit vector. The tetrad vectors e^ and 
their corresponding tetrad forms e° satisfy the orthonormality conditions 

e^.eb = ei'e.'^g^, = e^e^ ee e°^e V = 77"^ (1.1) 

where (7^1/ is the space-time metric, g^'' its inverse and 77"'' = r/ab the Minkowski 
metric. In the previous expressions and throughout this paper, we use Greek indices 
{fi, v, p,a . . .) for abstract indices and the beginning of the Latin alphabet (a, &, c, c? . . .) 
for tetrad labels. Since the tetrad labels run from to 3, we also use Latin indices 
starting from the letter i (that is z, j, k,l . . .) with values ranging from 1 to 3 to label 
the spacelike vectors or forms of a tetrad. We then reserve the label o for the timelike 
vector and form in a tetrad. A momentum can then be decomposed as 

p = p'^ea = p°eo p*e, , (1.2) 

where the components can be extracted as 

= p.e'' = eX ■ (1-3) 
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It can be decomposed into a magnitude p° and the a direction vector n according to 

=P°{e^ + n^), n.eo = n^e^ = 0, n^i^ = n.n = 1 . (1.4) 

This decomposition can be used to define a screen projector 

'5'eo^,.(p) = 9tii^ + e°^e° - , (1.5) 

which projects on a space which is both orthogonal to Gq and orthogonal to the 
direction n, since 

^e„^.e„^ = 0, 5e„^,p^=0, ^e„^X=0, 5e„ ; 5e„ - ■ (1-6) 

In this paper, a projected tensorial quantity is orthogonal to Bq only and a screen 
projected quantity is orthogonal to both n and So- If there is no risk of confusion about 
the choice of the observer in this screen projector, then we will use the notation 5^1, 
instead of Se^^y If we also omit the dependence of the screen projector in the photon 
momentum, then we abbreviate iS'^y(p) in 5*^1,, and S^^{'p') in S'^^. The polarization 
of a photon is represented by its polarization vector e(p) which is a complex spacelike 
unit vector (e*(p)e^(p) = 1) and taken in the Lorentz gauge (e^(p)p'^ = 0). Since 
there is a residual gauge freedom (of electromagnetism) in the choice of the polarization 
vector, we will work with the screen projected polarization vector S'^e''(p) which is 
unique. 



1.2. The (screen-projected) polarization tensor 

The radiation is represented by a Hcrmitian tensor valued distribution function (which 
is thus complex valued) , also called polarization tensor [531 IMl [Ml [S5] satisfying 

F^,{x^,p'') p^^F^,{x-\p-) = 0, (1.7) 

from which we can form the distribution function of photons in a given state of 
polarization e by 

F^,ix'\pne*^ip)e''{p). (1.8) 

Here the are the coordinates used to label points on the space-time manifold. 
Throughout this paper, indices which refer to these coordinates are A, B,C, . . . if they 
run from to 3. Furthermore, indices which are /, J, K, . . . run from 1 to 3, and the 
time component index is O. 

For a given electromagnetic plane wave with potential vector amplitude — yle^ and 
wave vector fc'^ in the geometric optics limit, this polarization tensor can be defined |66j 

by 

^^(P.p)^^^^^'^) = li^^fSU^ - p)A\,{k)et{k) , (1.9) 

where is the Dirac function of dimension n. This tensor F^^ should not be confused 
with the Faraday tensor. Since the remaining gauge freedom of electromagnetism also 
affects the polarization tensor, we can define the screen-projected distribution function 

by 

^,(x-4,p-) = SP^S:Fp4x'',p-) . (1.10) 

This tensor is no more dependent on the residual gauge freedom and encodes all 
the polarization properties of radiation as seen by an observer having a velocity e,,. 
Similarly to the definition of the screen projector, if the context requires it, we will 
use an index notation of the type f^^ to remind with which velocity, and thus with 
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which screen projector, it is defined. The screen projected tensor has four degrees of 
freedom which can be spht according to 

with e^i,cr = egep^i/cr, and where ep^i^a is the completely antisymmetric tensor. P^i^, 
which encodes the degree of linear polarization, is real symmetric and trace free, as 
well as orthogonal to and n, and thus encodes two degrees of freedom, which 
can be described by the Stokes functions Q and U [58]. I{x^,p°') and V{x^,p°-) 
are respectively the intensity (or distribution function) for both polarizations and 
the degree of circular polarization. We also define the normalized (screen-projected) 
polarization tensor by 

'^(^•^ - fa ~ J ■ y.^-^'^) 

1.3. Stress-energy tensor and energy-integrated functions 

For a distribution of photons, a stress-energy tensor can be defined by 

T^'^ix^) ^ eiie^ Ship.p)lix-\p^)py''^0^^ . (1.13) 

Performing the integral over p° (we choose the convention in which we integrate on 
the two mass hyperboloidcs), we eliminate the Dirac function and it leads to 

r-(.-4) ^ e^,-- (I /(,^p^)pV^^) , (1.14) 

where now the intensity distribution fimction has to be considered as a function of p', 
and p° is positive and taken on the mass shell, that is p° = \fp^ . and is thus identified 
with the energy of the photon. Splitting the integral over d'^p* into magnitude and 
angular direction leads to 

T^''(x^) (/ ^(^^pO {J>°f N''N''^P°A^^ , (1.15) 

where d^fi is the differential solid angle associated with the unit vector n% and where 
iV' = and iV° = 1. This motivates the definition of the (energy-) integrated 
counterparts of /, V and P^^ which are 

I(x^,p\n^Wfd.p\ (1.16) 
V(x-\p°yWfAp\ (1.17) 

P^.(x-4,p°,nO(p°)3dp°. (1.18) 

X is the brightness, V is the degree of linear polarization in units of I, and Vp,y is the 
tensor of linear polarization in units of X. 
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1.4-- Description of massive particles 

For massive particles such as electrons (e), protons (p) or cold dark matter (c), we 
do not need to describe polarization and thus we can rely solely on a the distribution 
functions gc, 5p and gc (chosen to describe the two helicities). Additionally, following 
common practice in cosmology, we will refer to electrons and protons together as 
baryons though electrons are leptons. This is motivated by the fact that most of the 
mass is carried by protons which are baryons, and the Compton interaction between 
protons and electrons makes these two components highly interdependent. For a 
particle with impulsion g", we use the following notation 

n' = 4' P^V^^, ^ = f3q\ n' = nyp, j = {l - p^'^^ (1.19) 

where now we have to distinguish between the unit vector and the velocity vector 
n' because of the mass m of the particles. The stress-energy tensor can be defined in 
a similar manner to equation (jl.isp by 

T'^-^ix^) EE ei^ei; (^J Shi^.ci m^)g{x-\q^)q^q''^0^^ . (1.20) 
Integrating over q°, this leads to 

T^^'^{x^) EE j^el^el^ (^J q')q"N-N'd^q'^ , (1.21) 
where we recall that N"" = (1, n'), and q" is taken on the mass shell {q° = ^ qiq^ + yn^). 



1.5. The fluid limit 



The stress-energy tensor of radiation or matter is equivalent to the one of an imperfect 
fluid with stress-energy tensor 



T, 



(1.22) 



In this decomposition p is the energy density, P is the pressure, is the fluid velocity 
and n^'', which satisfies 11^^ — w'^n^i,y — 0, is the anisotropic stress. For an isotropic 
distribution of radiation, that is where X(a;'^,p') depends only on the magnitude of p% 
it is straightforward to show by comparison of the expressions (|1.15p and (|1.22p that 
P = p/2i. Similarly, for a set of heavy particles (or non-relativistic particles), that is 
with \fcfqi <C TO, the pressure satisfies P p and the anisotropic stress tensor is also 
similarly small. For cold dark matter, it is assumed that the mass of particles is large 
enough so that we can use this approximation. 

In the case of electrons and protons, that is baryons, the Coulomb interaction ensures 
that the distribution of momenta follows a Fcrmi-Dirac distribution in the reference 
frame where they have no bulk velocity [2] , at least as long as the baryonic matter is 
ionized. This distribution is isotropic in this adapted frame and depends only on A 

s -1 



fffd(A) 



cxp 



^(A2-f m2)-^ 



T 



1 



(1.23) 
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As a consequence, the anisotropic stress vanishes, and in this adapted frame the 
baryons are then ideaUy described by the energy density and the pressure 

47r 



Pb - ^ j 5fd(A) v/A2 + m^X^dX , (1.24) 

If we can neglect the chemical potential fi, which is the case in the cosmological 
context, then for non-relativistic particles we obtain 

/ 3T \ 

1 + — Ub, (1.26) 
\ Zm J 

Pb^n^T. (1.27) 

The baryonic matter is ionized roughly until recombination where the temperature 
of photons is of order Tlss — 0.25 eV [S]. For electrons, the thermal correction is of 
order Tlss/^^c — 0.25/511000 ~ 0.5 10~^, and this ratio is even m^/me — 1836 times 
smaller for protons. We thus deduce that the baryons either have no anisotropic stress 
because of Coulomb interaction, or have a very tiny pressure and anisotropic stress 
after recombination has occurred. However, as discussed in section 16.21 this thermal 
correction is still of order of the metric perturbations and should not in principle be 
ignored for second order computations. Nevertheless, the thermal corrections are not 
relevant for computing the bispectrum and it is for this reason that we will, from now 
on, drop terms in T/m and describe baryons as cold matter (but not dark since it can 
interact with radiation) . It should be mentionned that for neutrinos these conclusions 
are not valid anymore since they are very light [67j . We will here assume that they are 
light enough to be treated as coUisionless radiation, and the equations which govern 
their evolution can be found by setting ctt = in the equations for photons, where 
ut is the Thomson cross section. 



1.6. Multipole expansion for radiation 

Functions of can be viewed as functions of (p°,n°) and we can separate the 
dependence into the energy and the direction of the momentum. The dependence 
in the direction can be further expanded in multipoles using projected symmetric 
trace- free (PSTF) tensors, where projected means that they are orthogonal to Gq. For 
instance, / can be expanded as 

oo 

I{x^y,n'^)=Y,IaJ.x''y)n^, (1.28) 

where the la,, = Iai...ai are PSTF. For the lowest multipole, i.e. the one corresponding 
to £ = 0, we use the notation /0. Note that we have defined the notation = 
. . .n"'. We also remind that n"^ = n.e° = n'^e'^. Since n is projected, n° = and 
thus if any of the indices in la,, is o, then the multipoles is chosen to vanish. These 
multipoles can be obtained by performing the following integrals 

where (. . .) means the symmetric trace-free part. A similar expansion can be 
performed on V, by replacing / by y in the above expressions, as well as for their 
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=T°°, (1.30) 

T' =47rAf^T°\ (1.31) 

= iTiA^^T^'^K (1.32) 

As for Pab = P^v^a^b 1 it can be expanded in electric and magnetic type components 
according to [SI EHl [Ml IM] 

Pab{x^,pl = £ [E,bc^[x^y)n^ - ny\^Bb)ac^[x^y)n^'\^ , (1.33) 

where the notation TT denotes the transverse (to n) symmetric trace-free part, which 
for a second rank tensor is 

[-''^afc]'^'^ = S^'^Sfj'^Xcd — -^SabS'^'^Xcd ■ (1-34) 
The electric and magnetic multipolcs can be obtained by performing the integrals 

E^ix^'y) = M^A^' J n^a_^Pa,_,a,){x'\p''y)d^n , (1.35) 



Ba^{x^,p°)=MfAj' J nbe'>\^na^P,,^_^^d{x-\p'',nnd'n , (1.36) 



where 



1.7. Transformation rules under a change of frame 

1.7.1. The photon momentum So far, everything was defined with respect to an 
observer having a velocity Gq. How does all this machinery transform when the 
radiation is observed by an observer with a different velocity Bq? Since two velocities 
can be related by a Lorcntz transformation, there exists a vector v such that 

eo = 7(e„ + v), 7= ^ = v.e° = . (1.38) 

V 1 — v.v 

We deduce immediately that the magnitude and the direction unit vector of the photon 
momentum transform as 

/ = p.e° = 7p° (1 - n.v) , (1.39) 

n = -(Go + n) - 7(eo + v) . (1.40) 

7(i — v.n j 

We remind that the direction, as observed by the transformed observer, is given by 
the decomposition p = p°{eo + n). These rules imply the following transformation 
rule for the screen projector 

27 / 7 \ ^ 

Sf_cu = Sf,,, + -^P(pSi,)pVP + f j P^PySafjV'^vf^ , (1.41) 

which implies directly the following useful relations 

S^y ~ S^Sj^ Spcr , n^^^pa = IT'^^papS pS^„ . (1.42) 



The radiative transfer at second order 



9 



The last relation can also be demonstrated easily by noting that n'^efj,af:iS°pS^^ is by 
construction orthogonal to Gq and n, and since it is also obviously antisymmetric in its 
two free indices, it has to be proportional to fi^i^pa- By contracting both expressions 
with themselves we obtain that they are indeed equal. 

The rest of the tetrad can be transformed without rotation, that is with a pure boost, 

by 

~e''^K\e\ ea^eXa, (1-43) 
and the components of this transformation are given by 

„2 



A° 



7 



(1.44) 



^ 1 + 7 " 

where we remind that i;' is the component of v along the tetrad e*, that is = v.e' . 
The transformation rule for the photon direction when expressed along tetrads is thus 



1 



7(1 — n.v) 



7 



(1 + 7) 



7W 



(1.45) 



1.7.2. The radiation multipoles It can be easily checked that for a vector orthogonal 
to p, such as the polarization vector, S^^S'^^e'^ = S''^e°'. As an immediate consequence, 
we deduce from equation p.9p that the screen-projected polarization tensor transforms 
according to 

/^.(a;-4,p°,n") = S^J''^Upix^,p'',n-). (1.46) 

We deduce from equation (|1.42p that transforms following the same rule, whereas 
I and V transform as scalars 

/(p°,n°) = /(p°,n°), f (/,n") y(p°,n"). (1.47) 

Here and in the rest of this paper, we omit the dependence in the position x-'^ to 
simplify the notation. We can deduce from equation (|1.39p that the differential solid 
angle transforms according to (this can also be deduced from using the transformation 
rule of p" and the fact that d'^p^/p° = p°dp°d^n is a scalar) 

1 2 



1 



dQ. 



7(1 — v.n) 

and this can be used to deduce the transformation rules of the multipoles 



(1.48) 



■{at) 



(1.49) 



/OO 
dn [7(1 - v.n)] -2 Kip'i^'i^ - "V)->^' 

In the previous integral, Ii,^, [p°7~^(l — n'^Vc)"^] has to be considered as a function of 
the direction n°. It is thus necessary to Taylor expand it as 

■ T 

where we define 



^ / — i C \ 

hAp'l-\l ~ n^v.)-'] = E A / V: ) ll:hp'l~'] , (1.50) 

— ^— ^ n! V 1 — 7 '-n'^Vr / — 

4;>(p°)^(p°)" (1.51) 



d{p°y 
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with the conventions /^^ = ij^p and I^^ = ij,^^ ■ Under this form, it is then possible 
to perform the integration using the following well known integrals [70] 
0, if fc = 2p+ 1. 

(1.52) 



47r 



S(i^^^ . . . S'^"-^-!"-^ , if k^2p. 
k + 1 I i 

Note that we have used the standard notation (. . .) for the symmetrization of indices 
which leaves unchanged an symmetric tensor and wc will also use the notation [. . .] for 
the antisymmetrization of indices which leaves unchanged an antisymmetric tensor. 
The integrals (|1.52p used in the transformation rule (|1.49p are ideally suited for a 
tensor calculus package, and we used xAct [7T] to compute them. 
Since the result of equation (|1.49p will involve terms like -f|7^[P°7^^] ~ Ibl^\p"{^ ~ 
n'^Vc)], it will require an additional Taylor expansion in order to have the result 
expressed only in function of the lj,'^^\p°] or lj,^\p°]- Note that the expression of 

iai(p°) in function of the ll,^^[p°], which is the choice that we make in the expressions 
that we report below, does not depend on the direction n', whereas it does when 
expressed in function of the I^^^ [p°] since p° is unambiguously defined from p° only 
once a direction is specified. A similar method, with similar definitions can be used 
to determine the transformation rules of the electric and magnetic multipoles [58] . At 
first order in the velocity v, we obtain the following transformation rules 

4 (/) = lajpl + ''^^f^^^'l^^Jp') + ^^/^'b^JP''^ 



E,,{P ) = %(P ) + (^+i)(2^ + 3) - E,J^P ) + JlTmM) '"'^ ^ 
-(£-l)«(,,i?,,_,)(/) + «(,,<,_^>(/) 



B (ri>\-P, . 5^ , (^-l)(^ + 2)(£ + 3) , , (£-l)(£ + 3) , , 
- (£ - l)i;<,,i?^)(/) + v^,fi'^_^{p') 

As for the energy-integrated counterparts, they transform at first order in v as 
and 

% - %- 1)^2^+3") "''^''f^-(^+3)^(a.ga^)+ (^^^^^%/a^)c .(1.58) 

We report in [Appendix A| the transformation rules up to second order in v for 
multipoles of further interest. 
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1.8. The Liouville equation 

In this section, we present the equation which governs the free-streaming of photons 
directly in the tetrad basis though nothing prevents it from being expressed with 
formal indices. The evolution of the polarization tensor is dictated by the Boltzmann 
equation 

L[fab{x'',p°,n')]^Cab{x-'',p",n'). (1.59) 

L is the Liouville operator whose action on TT tensors like the screen-projected 
polarization tensor fab is given by |58j 



L[fab{x^,p^)]^S:S,^ 



P ^hfcd{X ,p ) + 



(1.60) 



dp^ ds 

where s is the affinc parameter along the particle geodesic. V is the covariant derivative 
which in the tetrad basis is related to the partial derivative by 

p''VHfcdix'\p'')^p''dMx'\p'')+p''u;JfMix'\p^)+p'^u;Jf,b{x^,p^), (1.61) 

where the uiabc are the Ricci rotation coefficients (see |Appendix C| and [72] for details). 
Cab is the collision tensor whose expression will be detailed in section [S] 
In the case where the collision tensor can be ignored, that is when the collision of 
photons with electrons or protons can be neglected [the latter type of collision can be 
always ignored compared to the former since its cross section is reduced by a factor 
(too/ttt-p)^], this reduces to the Liouville equation. The Liouville equation arises from 
the fact that the polarization vector e of a photon is parallel transported and thus 
satisfies p^V ^e^ = 0, and then it follows from the construction p.9p of the (not screen 
projected) polarization tensor that 

/V.F..(..^p^) + ^^^l|lli^^=0. (1.62) 

By using the expression (|1.5p for S'-'" and the property F^^{x'^ ,p°-)p^^ = 0, we 
obtain directly from the previous equation that the screen-projected tensor satisfies 
Llfabix^jP"")] = 0. It can also be shown that the Liouville operator preserves the 
decomposition of /^^ in an antisymmetric part (V), a trace (/) and a symmetric 
traceless part {P^v), that is 

L[/,fc(a;■^p'^)] = ]^L[I{x^y)]Sab + L[Pab{x^ , p")] + ]^L[V [x"" y)]n' e,ab , (1.63) 
with the Liouville operator acting on a scalar valued function like / or y according to 

L[I{x'\p-)]^p'^d,I{x-\p'^) + ^^(0^^ . (1.64) 

To see this, we need only to use the property ujalbc] — ^abc and remark that 

SlS%p'^iu^,fSi = p^oJucfSlaSi^ = , (1.65) 

SfaS^^p^u^ncfef, (X Sf^^V^/e^, = e/[„e,f = , (1.66) 
where we have used the definition = n^ea^i,. 

Additionally, since the affine parameter s is a scalar, the transformation properties 
of the Liouville operator and the collision tensor under a local change of frame is the 
same as the transformation property of f^i, given in equation p.46p |58| . 
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2. Gauge transformations and gauge invariance for tensors 

The formalism presented in the previous section is very general and can be applied to 
the description of matter and radiation in any type of space-time. In the cosmological 
context it proves useful to use the high symmetries of the large scale space-time to find 
the solutions of the Einstein and Boltzmann equations using a perturbation scheme. 
We will review the standard perturbation theory for tensors in this section as well as for 
a scalar valued distribution function in section [3] and extend it to the TV distribution 
function. 

2.1. First- and second- order perturbations 

We assume that, at lowest order, the universe is well described by a Friedmann- 
Lemaitre space-time (FL) with Euclidian spatial sections. The most general form of 
the metric for an almost FL universe is 

ds2 = g^^dx^'dx'^ (2.1) 

= a{ri)^{ - (1 + 2$)dr;2 + 2uJidx^dTj + [(1 - 2*)(5/.7 + hij]dx^dx'^}, 

where 77 is the conformal time for which the corresponding index is O, and a{r]) is the 
scale factor. We perform a scalar-vector-tensor (SVT) decomposition as 

ojj = diB + Bi , (2.2) 

hi J = 2Hij + diEj + djEi + 2didjE, (2.3) 

where Bi, Ej and Hjj are transverse {d^ Ej — d^Bj ~ d^Hjj = 0), and Hjj is 
traceless {H^ j = 0). There are four scalar degrees of freedom (<I>, 5', B, E), four 
vector degrees of freedom {Bj, Ej) and two tensor degrees of freedom [Hjj). As 
we shall see in section 12.31 the perturbation variables live on the background space- 
time and thus their indices are lowered and raised by the background (conformally 
transformed) spatial metric and its inverse, that is with S^'^ and Sjj. Each of these 
perturbation variables can be split in first and second-order parts as 

W = W'-^^ ^W^^^ . (2.4) 

This expansion scheme will refer, as we shall see, to the way gauge transformations 
and gauge-invariant (GI) variables are defined. First-order variables are solutions 
of first-order equations which have been extensively studied (see [73] for a review). 
Second-order equations will involve purely second-order terms, e.g. VF^^^ and terms 
quadratic in the first-order variables, e.g. [M^^^-*]^. There will thus never be any 
ambiguity about the order of perturbation variables involved as long as the order of 
the equation considered is known. Consequently, we will often omit to specify the 
order superscript when there is no risk of confusion. 

As we shall see in section [^31 4 of the 10 metric perturbations are gauge degrees 
of freedom and the 6 remaining degrees of freedom reduce to 2 scalars, 2 vectors and 
2 tensors. The three types of perturbations decouple at first order and can thus be 
treated separately. As long as no vector source terms are present, which is generally the 
case when no magnetic field or topological defect is taken into account, the first order 
vector modes decay as a~^. Thus, we can safely discard them and set Ej^^ = Bj^^ = 0. 
In the following of this work, we shall not include first-order vector modes for the sake 
of clarity. We checked that our arguments and derivation can trivially (but at the 
expense of much lengthy expressions) take them into account. 
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In the fluid description, the four-velocity of each fluid is decomposed as 

u'^ = l{6-^ + V^). (2.5) 

The indices of V"^ are raised and lowered with the (conformally transformed) 
background metric, that is with rjAB and rj^^ . The perturbation V'^ has only 
three independent degrees of freedom since u must satisfy u^^u^ = —1. The spatial 
components can be decomposed as 

= d^V + V^ , (2.6) 

being the vector degree of freedom {diV^ = 0), and V the scalar degree of freedom. 
As for the energy density and pressure, similarly to any quantity which does not vanish 
in the background space-time, they are decomposed according to 

p = p + p^''> + ^p('^ + ..., P = P + P^^^ + ^P^^^ +... . (2.7) 

It is also common practice to define a conformally transformed anisotropic stress by 

T^CD^^^CD^ 7r'^ = n'^^, TTCD = a'ncD . (2.8) 

It should be mentionned that due to the symmetries of the background space-time, 
the anisotropic stress is already a perturbed quantity. At first order, the formalism 
developed by the seminal work of [7j provides a full set of gauge- invariant variables. 
Thanks to the general covariance of the equations at hand (Einstein equations, 
conservation equations, Boltzmann equation), it was shown that it was possible to 
get first-order equations involving only these gauge-invariant variables. In addition, 
if these gauge invariant variables reduce, in a particular gauge, to the perturbation 
variables that we use in this particular gauge, then the computation of the equation 
can be simplified. Actually, we only need to derive the equations in this particular 
gauge, as long as it is completely fixed, and then to promote by identification our 
perturbation variables to the gauge-invariant variables. Thus, provided we know this 
full set of gauge invariant variables, the apparent loss of generality by fixing the gauge 
in a calculation, is in fact just a way to simplify computations. Eventually we will 
reinterpret the equations as being satisfied by gauge invariant variables. The full set 
of first-order gauge- invariant variables is well known and is reviewed in [73j and |74j . 
As gauge transformations up to any order were developed, it remained uncertain |42j . 
whether or not a full set of gauge-invariant variables could be built for second and 
higher orders. This has been recently clarified [43], and the autosimilarity of the 
transformation rules for different orders can be used as a guide to build the gauge- 
invariant variables at any order. We present a summary of the ideas presented in |42j 
about gauge transformations and the construction of gauge-invariant variables |43j 
in a shorter version than in our previous paper [45] . A summary emphasizing the 
differences between the active and passive point of views can also be found in [75j . 



2.2. Points identification on manifolds 

When working with perturbations, we consider two manifolds: a background manifold, 
Aio, with associated metric g, which in our case is the FL space-time, and the physical 
space-time Mi with the metric g. Considering the variation of metric boils down to a 
comparison between tensor fields on distinct manifolds. Thus, in order to give a sense 
to "(;*^^^(P) = giP) — g{Py\ we need to identify the points P and P between these 
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two manifolds and also to set up a procedure for comparing tensors. This will also be 
necessary for the comparison of any tensor field. 

One solution to this problem |42| is to consider an embedding 4+1 dimensional 
manifold Af = A4 x [0,1], endowed with the trivial differential structure induced, 
and the projections Vx on submanifolds with Vq{N) = x {0} = A^o and 
ViiM) = M X {1} = Ml. The collection of A^a = Vx{M) is a foliation of TV, 
and each element is diffeomorphic to the physical space-time M. i and the background 
space-time A^o- The gauge choice on this stack of space-times is defined as a vector 
field X OYi Af which satisfies = 1 (the component along the space-time slicing R). 
A vector field defines integral curves that are always tangent to the vector field itself, 
hence inducing a one parameter group of diffeomorphisms ^(A, .), also noted 4'\{.), 
a flow, leading in our case from 0(0, p G Vo{N)) = p £ Vo{N) along the integral 
curves to 0(l,p G Voi-^)) = Q & Vi{JV). Due to the never vanishing last component 
of X, the integral curves will always be transverse to the stack of space-times and 
the points lying on the same integral curve, belonging to distinct space-times, will be 
identified. Additionally the property = 1 ensures that 4i\^x{'Po{N)) — V\{N)., i.e. 
the flow carries a space-time slice to another. This points identification is necessary 
when comparing tensors, but we already see that the arbitrariness in the choice of a 
gauge vector field X should not have physical meaning, and this is the well known 
gauge freedom. 



2.3. Gauge transformations and gauge invariance 

The induced transport, along the flow, of tensors living on the tangent bundle, is 
determined by the push-forward and the pull-back f72\ [75] associated with 
an element (jjx of the group of diffeomorphisms. These two functions encapsulate the 
transformation properties of the tangent and co-tangent spaces at each point and its 
image. Indeed, the pull-back can be linked to the local differential properties of the 
vector field embedded by the Lie derivatives along the vector field in a Taylor-like 
fashion (see [72] or [42]) 

k=0 

for any tensor T. The expansion of equation (|2.9p on VoiAf) provides a way to compare 
a tensor field on Vx{N) to the corresponding one on the background space- time Vq{N). 
The background value being Tq = C\T\-p^(^j^f^, wc obtain a natural definition for the 
tensor perturbation 

k—oo , ^ 

^^^^ - E -ki^'^^lw, = '^^•^(^^ ('-'^^ 

The subscript X reminds the gauge dependence. We can read the n-th order 
perturbation as 

^ y ,2.11) 

which is consistent with the expansion of perturbation variables of the physical metric 
in equation (|2.4p . since the physical space-time is labeled by A = 1. However, the fact 
that the intermediate space-time slices Vx (A/") are labeled by A removes the absolute 
meaning of order by order perturbations, as it can be seen from equation (|2.10p . The 
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entire structure embedded by f\f is more than just a convenient construction and this 
will have important consequences in gauge changes as we will now detail. 

If we consider two gauge choices X and Y, a gauge transformation from X to Y 
is defined as the diffcomorphism 

0x^y,A = {4>x,xr\4'Y.x), (2.12) 

and it induces a pull-back which carries the tensor AxTx, which is the perturbation 
in the gauge X, to AyTx, which is the perturbation in gauge Y. As demonstrated 
in [42l [62] this family (indexed by A) of gauge transformations fails to be a one 
parameter group due to the lack of the composition rule. It should be Taylor expanded 
using the so called knight-diffcormorphism along a sequence of vector fields ^n- For 
the two first orders, the expression of this knight-diffeomorphism is 

rv^T) = rx-^vAxiT) (2.13) 
= 0xa(T) + ^^^^<l^x.xiT) + ^{C^, + ClWxjT) + ... 

The vector fields fi, ^2 are related to the gauge vector fields X and F by ^1 = Y—X and 
^2 = [X, Y] . By substitution of the perturbation by its expression in equation (|2.10p , 
we identify order by order in A, and obtain the transformation rules for perturbations 
order by order. The first and second order transformation rules, on which we will 
focus our attention, are 

^T(2) - ;,t(2) = 2C^, xT(i) + (£5, + C'l )To . (2.14) 

General covariance, i.e. the fact that physics should not depend on a particular 
choice of coordinates is an incentive to work with gauge-invariant quantities. As we 
notice from equation (|2.14p . a tensor T is gauge- invariant up to n-th order if it satisfies 
xT^^^ = for any vector field ^ and any r < n, as can be deduced by recursion. A 
consequence of this strong condition is that a tensor is gauge-invariant up to order n 
if and only if Tq and all its perturbations of order lower than n either vanish, or are 
constant scalars, or are combinations of Kronecker deltas with constant coefficients. 
Einstein equation is of the form G — T = 0, and for this reason is totally gauge 
invariant. However, we cannot find non-trivial tensorial quantities (that is, different 
from G — T) gauge-invariant up to the order we intend to study perturbations, with 
which we could express the perturbed set of Einstein equations. 

Consequently, we will build, by combinations of perturbed tensorial quantities, 
gauge-invariant variables. These combinations will not be the perturbation of an 
underlying tensor. This method will prove to be very conclusive since a general 
procedure exists for perturbations around FL. Eventually we shall identify observables 
among these gauge-invariant variables and the fact that they are not the perturbation 
of a tensor will not matter. It has to be emphasized that the transformation rules of 
these combinations are not intrinsic and cannot be deduced directly from the knight- 
diffeomorphism since they are not tensorial quantities. Instead, we have to form the 
combination before and after the gauge change in order to deduce their transformation 
rules. 

We now summarize the standard way to build gauge-invariant variables up to 
second order. For simplicity we will not consider the vector part of the gauge 
transformations at first order, since we will not consider first order vector modes 
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in the metric and fluid perturbation variables (again, this could be done, but would 
just obfuscate the explanations). In the following, we split and ^2 

with djL^^"^^ = 0. We conventionally choose to lower and raise the indices of d^L + 
with the (conformally transformed) background metric, which implies that 

6/ = a'diL'-^^ , 6/ = a^diL'-'^ + a^lf . (2.16) 
2.4^. First- order gauge-invariant variables 

In the subsequent work we present the transformation rules of perturbed quantities 
in a simplified notation. Instead of writing yW^^^ = jfVF^''-' + / (^1, ..,^r), in order to 
state that the difference between the expression of the r-th order pertubed variable W 
in gauge Y and in gauge X is a function / of the knight-diffeomorphism fields ^1 , 
we prefer to write W^"^^ VF^''-' + / (^1, .., £,r)- We remind that the expressions of the 
fields (Cri)i<n<r necessary for the knight-diffeomorphism are expressed in function 
of the gauge fields X and Y [see below equation (|2.14p ]. From the transformation 
rules (|2.14p we deduce that the first-order perturbations of the metric tensor (|2.ip 
transform as 

$(1) ^ $(1) +r'(i) + 7^r(i) (2.17) 

^ ij(i) _ T(i) L'W (2.18) 

^ ^r(i) -T^T^i) (2.19) 

^(1) ^ £;(!) (2.20) 

ijjy ^ h\'} , (2.21) 
while the quantities related to matter transform as 

p(i) ^ + 

V'-^^ -^V^'^-L'^'^ (2.22) 
^IJH) ^ ^/./(i) , (2.23) 

where a prime denotes a derivative w.r.t. to conformal time 77, and where Ti. = a' /a. 
From now on, we shall refer to these first-order transformation rules defined by 
as T^, ($(!)), T5,(B(i)),... or simply T($(i)), T(B(i)), ... For instance r($(i)) = 

$(1) + nT''^\ 

We first note that the first-order tensorial modes and the first-order anisotropic 
stress are automatically gauge invariant. We can define gauge invariant variables, 
by transforming and 5'^^-' towards the Newtonian gauge (NG) [10]. This 

transformation is defined by the vector field C^nq decomposed in T^^nq = B^^'^ — 
E'^^\ -^'^NG = ^E^^\ and it transforms the perturbation variables as 
B^i) (2.24) 
E^^^ (2.25) 
^(1) ^ 1,(1) = ^^<j,(i) ^ $(1) +7^ - + - iJ'd))' (2.26) 

^(1) ^ ^(1) = ^^^(1) ^ ^(1) _ ^ (^^(1) _ ^/(i)^ _ (2.27) 
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Similarly the gauge-invariant variables that would reduce to 6p, 5P and v are 

^(DP = ,^P(1) - P(i) + P' f P(i) - 



7r"(i) = ^^7r"(i) =7r^''(i). (2.28) 

Since we have ignored the vector gauge degrees of freedom, B^^^ and E^^^ are 
the two gauge variant variables of the metric perturbation while l>'-'^' and 'I''^^ are 
the gauge-invariant part. As mentionned before, we then force the gauge-invariant 
variables in the perturbed metric by replacing with "l'^^ — H, (P'^^ — P''^^) -f 

(P^^^ — P'*^^^) and applying similar procedures for p^^'> and P^^\ When 

developping Einstein equations, we know that general covariance will eventually keep 
only gauge-invariant terms. Thus, we can either do a full calculation and witness 
the terms involving the degrees of freedom P^^^ and P^^^ disappear, or perform the 
calculations with P^^^ and P^^^ set to zero and obtain the perturbed Einstein equations 
only in function of gauge-invariant variables. The latter simplifies the computation, 
which is useful when going to higher orders. This procedure means that we decompose 
the perturbed metric in a gauge-invariant part and a gauge variant part as 

g(i)^g(i)+£ 5, (2.29) 
as it can be seen from the transformation rules under a gauge change characterised 

by 6 

-f^NG^ -C^lc+ei, (2.30) 

and that eventually, only g^^'^ will appear in the equations. This property which is not 
general but happens to hold in the case of cosmological perturbation (i.e. around FL 
metric) is the key to extend this construction to second order. 

It should be noted that this procedure, although achieved by defining gauge 
invariant variables which reduce to the perturbation variables in the Newtonian 
gauge, can be extended to other types of gauge-invariant variables which reduce to 
perturbation variables in another gauge. 
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2.5. Second- order gauge-invariant variables 

For second-order perturbations, equation (|2.14p gives the following transformation 
rules 

$(2) ^ $(2) + y'(2) _^ -^y(2) ^ 
5(2) ^ ^(2) _ y(2) ^ ^/(2) ^ 
4r(2) ^ ^(2) _ ^7^(2) ^ 
^(2) ^ i5(2) + ^(2) + 

H-(2) ^ 7:r(2) , c, 

^ ^ -,y(2) ^ 5^ 

p(2) _^ p(2) _^ p'p(2) ^ 
y(2) _ ^(2) _ p'(2) + 
^/(2) ^ yl(2) _ ^,/(2) ^ 

^"(2) ^^/./(2)+5^„^ (2.31) 

where the source terms are quadratic in the first-order gauge transformation variables 
T(i),L(i), and the first order metric perturbations <^^^\<^'-'^\B'^'^\E^^'> and E^^] . We 
collect the expressions of these terms in [Appendix B[ In the rest of this paper, we 
shall refer to these second-order transformation rules associated with (^) = (^1,^2) as 
7^j)($(2))^7-^^^(5(2))^ or simply r($(2)) 7-(p(2))^ . These transformation rules are 

much more complicated than their first-order counterparts. However, the combination 
defined by F = (7^2) _|_ 2Cai) g^^'^ -\- C^n) g enjoys the simple transformation rule 

«^NG 4_NG 

F F + L , r,(i) , ,g under a gauge change defined by ^2 and ^1 (see [15]). As 
a result, its transformation rule mimics the one of first-order pertubations under a 
gauge change. This means that if we decompose F in the same way as we did for the 
metric with 

Bf ^B(2)+5B(eli^NG) 
Ef ^ E'^^^ + SEiitio) 
Bf, = Bf^ + Sbi (C*^ng) 
Ef, = Ef^ + Se, (C^ng ) 

Hp, J = Hj^j^ 4- Sh,j{^'lJmg)^ (2.32) 
then the transformation rules for these quantities will be similar to those of 
equation (|2.24|) . but with the vector ^2 + K-»nGi^i] instead of ^1. Consequently, we 
shall use the same combinations (taking into account the vector contribution at second 
order since it is does not vanish at this order) in order to construct gauge-invariant 
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variables which are 

$(2) =<S>p + (Bp - E'p)' + n{BF- E'p) 
^(2) ^-^p-niBp- E'p) 
^ Bp, - E'p, 

Hf]^Hp,,. (2.33) 

This procedure is equivalent to transforming quantities in the Newtonian gauge since 
it transforms B, E and E^ into a null value up to second order. This transformation 
is defined by S,\,^g ^^^^ decompose in 

Ti^) ^ ^ B(^) - E^^ + Ss (eaio) - Sp 



+NG — ^ ^E (^S^NG 

l'SL = - E^''' Sp: (e^) . (2.34) 

The second-order gauge-invariant variables can thus also be defined by 
$(2) = $(2) 

^ — NG ^ 



^(2) = ^(2) 

^ — NG ^ 



$(2) ^ ^(2) 

— NG-^I 

Hf] ^ ,,Hf] (2.35) 

= n(2) 
p(2) ^ p(2) 

V^(2) = ^(2) 
— NG " 

- -f{2) 

^"(2)^;^^/.7(2), (2.36) 

where the index NG indicates that wc transformed the quantity with the formula 
(|2.14p . with the vector fields ■f^^NQ and ■C^nq defined above. This procedure means 
that we have split the second-order metric according to 

g(2)=g(2)+£ (2, g + 2C a, g^^^ - g, (2.37) 

"S — NG "S^NG ~^^NG 

where g^^'^ is the gauge-invariant part and — 'CI^nq the gauge variant part, as it can be 
seen from the transformation rules under a gauge change characterised by (^1,^2) 

^(2) 

-ei!l,G ^ - ^^'ng + 6 + [e!!i,G> Ci]- (2.38) 



3. Gauge transformations and gauge invariance for a distribution function 

3.1. pre-Riemannian distribution function in the unpolarized case 

So far, we have set up the mathematical framework to identify points between the 
background space-time and the perturbed space-times through a gauge field X. This 
enabled us to define the perturbation of tensors and to calculate their transformation 
properties under a gauge transformation. However this allows only to perform a fluid 
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treatment of the radiation. In the statistical description for a set of particles, we 
assume that each particle has a given impulsion and is located at a given position 
[75] . The equations then have to describe the phase space distribution of the particles. 
If the number of particles is high enough, we can define a probability density, the 
one-particle distribution function, of finding a particle in an infinitesimal volume of 
the phase space. Now, let us focus our attention on this distribution function. The 
distribution function is a function of the point considered (i.e. its coordinates x"^), 
and also a function of the tangent space at this point whose coordinate we label by 
p'^du- There is no special reason for this function to be linear in p'^d^, but we can 
expand it, without any loss of generality, in power series of tensors according to 

/(x^,p'')=5]F^,..^,(x-'^K\..p^^ (3.1) 

k 

Here / stands for either the intensity /, the degree of circular polarization V or the 
matter distribution functions g^^ gb and g^. From the previous section we know the 
transformation rules for the tensorial quantities -F^i..^;,, thus / transforms according 
to 

% [fi^^^P")] ^E^(C) [F,...,A^^)]p'''-r\ (3.2) 

k 

where 7(^) refers to the knight-diflFeomorphism with the set of vectors (^1,^2, •■•)■ 

As we do not necessarily want to refer explicitly to the decomposition in 
multipoles, we use the fact that for any vector ^ = ^^d^, which defines a flow on the 
background space-time Vo{Af), we can define an induced fiow (a natural lift) on the 

vector tangent bundle TVoiM) directed by the vector field = ^^df^p^id^^f")-^ . 
This implies the useful property 

(^^^...Mjp^^-y'" = ^T? (F^^.^.p^^.p^") . (3.3) 

With this definition, we can rewrite the transformation rule for / as 

='r^r,)[f{x-\p'')], (3.4) 

where now T(t^) refers to the knight-diffeomorphism with the set of vectors 
(T^i, T^2j ■•■)■ The evolution of the distribution function is dictated by the Boltzmann 
equation, and its collision term can be expressed in the local Minkowskian frame 
defined by a tetrad fields e^, from known particles physics. For this reason, the 
framework developed to define gauge transformations for a general manifold has to be 
extended to the case of Riemannian manifold. Instead of using the coordinates basis 
Oa to express a vector of tangent space as p = p'^Oa, we prefer to use the tetrads 
basis Ba and write p ~ p^'Sa- In terms of coordinates, this means that the distribution 
function is a function of x"^ and p", and that we favor the expansion (|1.28p instead 
of the expansion (|3.ip . When expressing the physics with the tetrad fields, the metric 
is not just one of the many tensors of the theory whose properties under a gauge 
transformation we need to know, but rather a central feature of the manifold, since 
it determinates the tetrads (up to a Lorentz tranformation) required to express the 
distribution function. 
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3. 2. Perturbation of tetrads 

3.2.1. Pulling hack the tetrad With the formalism developed for tensors, we carry 
this tetrad field onto the background space-time using a gauge field X with 

k—oc ^f. 
k=Q 

xe^-^^C^e e.EE^ei"), (3.5) 

and similar formulas for e°. 

As is a basis of the tangent space on the background space-time (and a 
basis of its dual space), x^a and xe" which also lie on the background space-time can 
be expressed in the generic form 

X^a = xRa^b , X^'' = 6° xSj* , xRa xS^ = xS^ xRq = (^ai (3-6) 

where 

xRab = X! fc[ ^^^ab ' xSab = "fcT ^^^"^b ' (3-7) 

fe ■ fe ■ 

Order by order, this reads 

xei") = xi?i"^'efc, A-e'(") = (3.8) 

3.2.2. Normalization condition Tetrads are four vector fields which satisfy 
equation (jl.H) and are thus related to the metric. Consequently, the perturbations of 
the tetrad defined above are partly related to the perturbations of the metric. When 
pulled back to the background space-time, equation (jl.ip implies 

4>xx{'nab) = Vab = 0A,x(e^^e^5^^) 

= <l>lxKWx.x{e''bWx.x{9^^)- (3.9) 

Identifying order by order we get in particular for the first and second orders 

Gfc.A-ei'^ +e,.Aef' + ^.g^^^ (e,,, e^) + ^ei'^ xe^ (3.10) 
+ J'^ [xei'\e,) + J'^ (e„,Aei^))=0, 
where a dot product here stands for g { , ). From the constraints (j3.10p , we can 

(n) 

determine the symmetric part of i?),^ as 

xR'il)= -lj'Hea,e,), (3.11) 

« J " x..a X..,, , 
which are related to the components of the inverse by 

x^il^ = - xRH^ , (3.13) 
x^lf = - xRZ^ + 2xR'^^^xRil^. (3.14) 



5'^' (e„ xRi'Je') - xR^^^'xRi'J , (3.12) 
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The antisymmetric part, xRiab], still remains to be chosen as it corresponds to 
the Lorentz transformation freedom (boost and rotation), which is allowed by the 
definition (jl.ip . A first and easy choice would be xR^^b] ~ ^ ^"^y ^- However, as 
mentioned above, we eventually want to decompose a vector p^dA on tangent space 
as 

p''dB=v''ea=p''ej^dB. (3.15) 

and identify p° with the energy and with the momentum (although conserved 
quantities are generally ill-defined in general relativity, energy and momentum can be 
defined when performing perturbations around a maximally symmetric background 
[77] as it is the case here). When working with coordinates, we want to express 
physical quantities, as measured by a set of fundamental observers that we need to 
choose. If we choose the velocity of these fundamental observers as orthonormal to 
the constant time coordinate hypersurfaces, then we require |78j e° ~ dry, which is 
equivalent to choose xS*'"' = for any n. This choice allows us to fix the boost in 
by imposing the condition xS\^)i^ = — xS\^l^ = — xS^^}-^ , and it can be checked by 
recursion that this implies 

xR^ = 'xI^=-xI^:ly (3.16) 

We also fix the rotation by requiring xS'^j"] = 0, and it can be checked similarly, that 
this imphes = 0. 



3. 3. Gauge transformation of tetrads 

Under a gauge transformation, we can deduce the transformation properties of the 
tetrad from those of the perturbed metric. In the FL case, we use a natural background 
tetrad associated to Cartesian coordinates 

eo = {do)/a, ei = {di)/a, (3.17) 

in order to evaluate equation p. lip . We report the detailed expressions for the 
transformation of the tetrads for the first and second orders in [Appendix C[ Since 
the tetrads are constrained by equation (|3.9p , the transformation is necessarily of the 
form 

r(e,) = Aj'rk„ight(e6) , Tie'') = Af,"rk„ight(e'') , (3.18) 

where Tonight refers to the transformation that one would obtain by considering each of 
the four vectors fields in the tetrad as independent vector fields, whose transformation 
rule is calculated using the knight diffeomorphism of equations (|2.14p . As a result 
of the choice = S'[.y] = 0, we get that A/ = , and the three vectors in a 
tetrad transform as independent vector fields. Similarly, if a vector field W satisfies 
W.Go = 0, then transform as a scalar field under gauge transformations (recall that 
= e'.W). This property extends of course to higher rank tensors, which means 
that the tetrad components of a projected tensor are to be considered as scalars under 
gauge transformations. 
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3.4- Gauge transformation of a distribution function 

3.4-.1. The perturbed distribution function with a tetrad Now that the transformation 
properties of the tetrads are known, we turn to the general transformation of a 
distribution function /(x'^,p°). First, by using the tetrad we can deduce a muhipole 
decomposition of the type p.28p from the decomposition (|3.ip by 



^f{x^,pn- (3.19) 
Since the momentum is constrained to be on the mass sheU, only three of the four 
components are independent, and we can thus choose the multipoles Fai{x^) to be 
projected, that is to say if any of their indices is o, then they vanish. Note that in the 
case of radiation, that is for f = I, then Fa,, = {p°Ylai ■ 

We do not need any additional identification procedure for the tangent spaces through 
a gauge field, in order to identify points of the tangent space of the slices TV\{N). 
Indeed, once the metric and a gauge field X are chosen, there exists a natural 
identification with the tetrad fields. First, and as mentioned before, we identify 
the points of M which lie on the same integral curves of X, that is, we identify a 
point P € VoiN) and ^x.x{P) € 'P\{M). Then, we identify vectors of their respective 
tangent spaces, if the coordinates of these vectors in their respective local tetrad frames 
Ga and Gai are the same. To be short, we identify p'^Ga and p^e^. The multipoles 
Pai{x^) are scalar fields, and they can be pulled back on the background space-time 
using the gauge field X, and we define in this way perturbations 

F,,{x^) ^ xFaAx-^) ^ V - xFif{x^) . (3.20) 



A 



This perturbation scheme induces a perturbation procedure for the distribution 
function / as 

A" 

xf{x^,P^) ^ ^^x/(")(x-4,P°), 



X. 



■f^"\x^,p^) EE ^ ^^F^){x^)p^. (3.21) 

e 



It is essential to stress that is not a perturbed quantity, it is a coordinate of the 
locally Minkowskian tangent space. However, the tetrad field allows us to see p as 
a perturbed vector since p(p") = Gap". In other words, for a given there is an 
associated vector whose order by order perturbation in a given gauge X is given by 



3.4-2. Gauge transformation of a distribution function We can deduce the 
transformation rule under a gauge change directly on the form p.lOp . pulled back 
to the background space-time, 

r [xfix^\pn] ^Y.'^[x:f,,...,,{x^)] r (xO . . .r (xe^;)^"^ . . .p"^ (3.22) 
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The first factor in this expression is tensorial. Exactly as for the pre-Riemannian case, 
its transformation rule is dictated by the knight-diffeomorphism, whereas we get the 
transformation rules of the tetrads from equations (jC.2p and equations (|C.4[) . As we 
do not necessarily want to refer explicitly to the multipole expansion, the first factor 
is rewritten by considering / as a function of p'^ using p° = x^Zp'^j applying 
equation (|3.4p . We then have to consider the resulting distribution function as a 
function of p", knowing that the inversion is now given by p{p"') = T{ea)p°'. This will 
account for T ( ^e^^) . . . T ( x^''^") in equation (|3.22p . In a compact form it reads 



T[jix-\p'^)] =7^^,){,/[x^eX]} ^. (3.23) 

To obtain an order by order formula, we explicit these three steps using a Taylor 
expansion. First, we use that 



exp(e>^x5,° 



d 
dp" 



{x'\ey)^ ,g{x'\pn, 



(3.24) 



in order to consider / as a function of p^. We then Taylor expand again the result of 
the knight-diffeomorphism in order to read the result as a function of p° , 

d 



T[Jix-\p'^) 



exp(etp"r(xi?.'')^)r(^0 (.5) 



(^-^e-^^p"). 



(3.25) 



The derivatives in the previous expressions have to be ordered on the right in each 
term of the expansion in power series of the exponential. When identifying order by 
order, we need to take into account the expansion in Rat and Sab, in the exponentials 
and also in the knight-diffeomorphism. The detail of the transformation properties 
of the first and second-order distribution function can be found in [45| . We report 
here the final result of the transformation rule for the first order and the second order 
distribution function for radiation 



T 



(1) 



- f(i) 



—p n diT + T— 
op° or] 



(3.26) 



T 



.(2) , ^(rp(2) 



TT' + diTd'L) 



^ {p*a/r(2) - 2p^ [{didjE + Eij + didjL) d^T - ^djT] 



-p°diTd'T + p\TdiT)' +p'di {d-'LdjT) + 2^p'diT} 



d{p°) dridp° 



Tp'diT + —4tT 



2^ — p^djT 



drj^ 



2^^^p°d^T -f 2d'LdiJ^^^ 



2T 



(3.27) 



dp" ' dp^ ^ drj 

In order to derive this formula, no use has been made of the normalization of the 
momentum. Consequently in the above transformation rule, / stands either for /, V 
or for ge, gh (where we conventionally prefer to use the notation q instead of p for 
the momentum of particles). Note here that the Einstein implicit summation rule still 
applies on index of tetrad type {i,j,k...) when contracted with indices of coordinate 
type (/, J, K . . .). This type of contraction arises from the fact that we have chosen a 
background tetrad field adapted to the coordinate system, that is with ^ and 
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3.4-3. Gauge transformation of the linear polarization tensor Now if we want to 
describe polarized radiation, it is straightforward to generalize this calculation. Indeed 
it is projected, that is foi = fio = foo = 0, and thanks to the remark at the end 
of section 13.31 each of the nine components /y can be treated like a scalar valued 
distribution function, and the same property is valid for P^. Furthermore, since 
for symmetry reasons the radiation is not polarized at the background level, that is 
Pij =0, its transformation rule reduces to 

- x^ii ' , (3.28) 



^ [.P^^) = .pS' + ^^^P'djT (3.29) 
a p(i) p(i) 

We can then form first-order and second-order gauge invariant intensity function, 
degree of circular polarization and tensor of linear polarization, by transforming these 
quantities in the Newtonian gauge 

/« EB,jWEBr^a,^Jx/(i)) (3.30) 
yW^,,FW^r_a, (xV^'A (3.31) 

P^P ^ .oPj ^ ^ju)^^ ( xP^^) , (3.32) 

(3.33) 



/(^) ^,,/(2)^r.., JxI^'A (3.34) 

y(^) ^ ^ r. JxV^'A (3.35) 

1 C^NG' ? — NG \ / 



p(2) ^ p(2) ^ ^ 

^t] — NG-^ij — (^f^'j, 



NG- 



,^^{xPS')- (3.36) 



The energy integrated counterparts '))'■"•' and ■p|"'' with n = 1,2, follow from a 
definition similar to equations (|1.16p . Since the transformation rule (|3.27p is valid also 
for the electrons and protons (or for baryons which is the sum of these two components) 
distribution functions, the corresponding gauge invariant distribution functions can be 
defined in a similar manner and we name them g^^^^ and g^^^^ for the first and second 
order respectively. 



4. The gauge invariant Liouville equation for radiation 

Now that we have a complete formalism to handle the perturbations of the fluid quan- 
tities and of the distribution functions, we will use it to express the perturbations of 
the dynamical equations satisfied by these quantities. In this section we focus on the 
perturbative expansion of the Liouville operator and in section [S] we detail the per- 
turbative expansion of the collision term. In both cases we will recover the standard 
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first order result and the novelty lies in the second order perturbation. 

In practice, we want to express the Boltzmann (or the Liouville equation if 
relevant) in function of rj rather than the affine parameter s, since we want to perform 
an integration on coordinates. We should thus multiply equation (|1.59[) by ^ = l/p'^- 
However, there is no point multiplying with the full expression of ^ since when 

, (fc) 

< k < n will multiply the {n — fc)-th order of the perturbed Boltzmann equation 
and these terms will thus vanish. This is the reason why we will instead multiply only 



focusing on the n-th order of the Boltzmann equation, the perturbations ( 4^ ) with 



by ( 



— (l/p'^"')^"'' — a/p", and we will use the notation 

(0) ,. 



, 0*B.C..(^) . (4.1) 

In order to compute the Liouville operator, we need to calculate the perturbed 
expression of and . They can be obtained from the geodesic equation 

^+c.6,,py = 0, with ^=/56p,=A6^ac(Pa), (4.2) 
as as 

that we pull back to the background space-time in order to extract order by order 

equations 

^) n = 0,l,2, (4.3) 

and from the pcrturbative expansion of 

^) n=.0,l,2. (4.4) 



It should be noted that the transformation properties of [] under a local change of 
frame is not similar to that of due to the factor l/p° which transforms according 
to equation (|1.39p [5B j. 

The result of the Liouville operator up to second order is of the form 

L[X,X(i),X(2)] =Z[X]+i(i)[X,X(i)] + i (l(2)[^^^(2)]^^(i)(i)[^^^(i)]^ ^ (4 5) 

where X stands for cither /, V or Pab- Similarly to Einstein equations, the second order 
Boltzmann equation will involve either terms which are linear in purely second order 
quantities that we gather in or terms quadratic in first order quantities that we 

gather in The linear dependence in these purely second order quantities is, 

as usual, the same as for the linear dependence obtained for the first order equation 
with respect to the first order quantities, and this means L'-^-'O = i'^-'O- In principle, 
once presented the expression of we only need to report the expression 

of L^^'>^-^'>[X , X^-^")] to fully express the second order equations. In practice however 
we will report both since we chose to neglect the first order vector degrees of freedom 
but we cannot neglect them at second order and they will thus contribute to -^^^-'O- 
Additionally, we also choose to neglect the first order tensorial perturbations. Though 
they are generated in standard models of inflation, their amplitude is expected to be 
much smaller than first order scalar perturbations, and furthermore they decay once 
they enter the Hubble radius. Similarly to vector modes, we will not neglect the tensor 
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modes at second order and their contribution will be also reported in i'^^O- 
Since the expressions in the decomposition (|4.5|) have a dependence in n'), we will 
also perform a multipolar expansion according to equation p.28p for the equations on 
/ and V (with multipoles L[I]ai and L[T^]a^) and according to equation p.33p for the 
equation involving Pat (with electric and magnetic type multipoles being respectively 
L[E]a^ and 

We will also define the energy integrated Liouville operator by 

C[]^j^l mp-fdp", (4.6) 

and use similar definitions for L'^l], Cab and C^, and this will lead to define multipoles 
for these quantities, with an obvious choice of notation. We recall that since the 
polarization tensor is a function of the position coordinates , that we do not write 
explicitly, and of the momentum components p", then the Liouville operator and the 
collision tensor have a dependence in the same variables. However, in order to simplify 
the notation, we will not write the dependence in p° throughout this section, which 
is dedicated to the Liouville operator for radiation, and throughout section [51 which 
is dedicated to the Liouville operator for baryons, but we will restore it in section [6] 
dedicated to the collision tensor. 



4-1. The background Liouville operator 

At the background level, space is homogeneous and isotropic. Consequently, the 
radiation is not polarized and its distribution function depends neither on the direction 
of the photon nor on the position in space . It only depends on p° and 77, which 

dl ^ dl 

dn'- djr 



implies that -M^ = = Pat = 0. Since the background geodesic equation implies 



(0) 

— —T-Cp", the Liouville operator reads at the background level 

where we recall that the derivative w.r.t. 77 is to be taken at p° fixed. Its energy 
integrated counterpart reads 

c*[T] = ^ + Ani. (4.8) 

orj 



4-2. The Liouville operator at first order 



The equation required to express the Liouville operator is the first order of equation l4.3 
which in the Newtonian gauge leads to 



dp° 
d?7 



(1) 



= P 



(4.9) 



and we obtain that the first order perturbed Liouville operator can be expressed in 
function of gauge invariant quantities by 



Wji^^'> - Hp° 



dj] ' dp° 

Its energy integrated counterpart is thus given by 



£#(i)rt jd) 



drj 



(4.11) 
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It will prove more convenient to decompose it in multipoles which are given by 



mx 



(1) 



:d''^)t:+diiji:>,^ (4.12) 



(2^ + 3) 

As for V, it follows the same equation as / but we will see in section [6] that it is not 
excited by the collisions and remains null. 

Following the same method for the linear polarization tensor leads to 



dP' 



(1) 



dr] 



+ n^djP, 



Hp 



QpW 

dp" 



dV. 



(1) 



(1) 



(4.13) 



(4.14) 



We can extract the electric and magnetic type multipoles of the energy integrated 
Liouville operator, and they read 



£#(l)[^(l)] 



(1) , ^^c(l) I (^-l)(^ + 3) ;,(!) 



(4.15) 



(^+1) 



£#(i)r^(i)i . _ , 4-^^(1) , (^-l)(^ + 3) ^(1) 



(^+1) 



(curl£(i)),,, 



where 



(curlX)^^, 



(4.16) 



(4.17) 



However, it can be shown that when the first order tensor and vector modes are 
neglected, which is our case in this paper, the first order magnetic type multipoles are 
not excited and remain null [15]. We will thus discard them from the computation 
in the following of this paper in order to simplify the intricate expressions of second 
order perturbations. In order to be consistent with equation (|4.16p we will also have 
to neglect (curlf^^))^^. 



4-3. The Liouville operator at second order 



Following the same method as for the first order case, we obtain the second order 
evolution of the energy p° which in the Newtonian gauge reads 

[-71^9/$(2) + Vl/(2)' + [djBf - Hf]' 



dp" 
drj 



= P 



+2($ - 5')n*a/$ + 4*^''] . (4.18) 
We also need the first order photon trajectory which in the Newtonian gauge reads 



i\ (1) 



(1) 



(4.19) 
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and the first order Icnsing equation whieh gives the evolution of the direction 



dn 



(1) 



dr] ) 

We then obtain the gauge invariant second order Liouville operator 



drj 



dp" 



(4.20) 
(4.21) 



^i 5^1,(2) + ^(2)' _^ (di^'^j^ - H 



r(2)' 
'iJ 



,9/ 

dp° 



n'di^ + ii' + 2($ - 'I')7i*9/$ + 4**' 



-2<|.L#(i)[/,/(i)] 
Its energy integrated version is 



dri'- 



,dl 

dp° 



drj 
-41 



„J9,$(2)+^(2)' + /'a^$(;)_^(2)' 



(4.22) 
(4.23) 



i ? 



4X 



2($ - *)n*a7$ + 4#*' 



-2$£#(i)[X,f(i)]. (4.24) 

We recall again that we have omitted some exponents stating that some quantities 
are first order perturbations variables, since at second order, the terms involved are 
either involving purely second order perturbations or terms quadratic in first order 
quantities, and there is thus never any ambiguity. Note also that, if L^[] and had 
been defined with (ds/drj) instead of (ds/dr/)'"' in the definition (|4.1[) . the last line of 
equations (|4.22p and (|4.24p would not be there, and this enables us to compare with 
the results of [IHl US] where such different choice has been made (though it is then 
inconsistent with the collision term of [JS] which has been calculated with (ds/dr/)'"' , 
but not with the collision term reported in [49j which seems to be defined with (ds/d77) 
though it is stated that it is the same expression as in [15j). We can then extract the 
multipoles of these equations and we obtain 



r'(2) 



mi. 



(2) 



./-f(2) 



-d'X 



li ' (2^ + 3) " ^ 



(4.25) 



4<5|J 



(2) 



H 



(2)' 
I1I2 



(2^ + 3) 



- 2{£ + 2)9-'($ + f) + 2^'a'^] X 



(7,<i> + 2i£ - 1)9(7, ($ + «-) + 2*9(7,] 2]^) 
45jl [2{^ - 2$)97$] + 4(5^X [(2$ - 4^)^'] . 



(4.26) 
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Following the same method as for the intensity part of the Liouvillc operator, we 
obtain the Liouville operator for the linear polarization tensor 



dP, 



(2) Qf)(2) 



i#(l)(l)[p(l)]^^a^^ 



dP. 



(1) 



ah 



and its energy integrated counterpart 



#(2) [73(2)1 



cd 



Qa Qb 
>J(7 



dV 



(2) 
ab 



drj 



^(2) 



2(l> + *)n'97p^p - 2S'^ (dj^ + 9,/$ 



(1) 



dri^ 



'9/$ + *'] -2$/:#(i)[75^^'] 
We extract the electric and magnetic type multipoles 



(4.27) 



(4.28) 



(4.29) 



(4.30) 



2 



(2) 

' (£+l)(2^ + 3)" "^^^ ' "'•''^^i^) 



(curlS(2))^ 



Sa '-I - 2(£ + 2)d \^ + #) + 2#a'^ 



(2^ + 3)(^+l) 
+ [8(9(j,$ + 2(^- l)9(j^,(l' + *) + 2#a(j,, 

- 8^'£„. - 2$f - snM,, , 



(4.31) 
(4.32) 



£#(2) [^(2)]^ 



(£+ l)(2£ + 3) 



(curl £(2))^ 



(4.33) 
(4.34) 
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4-3.1. Dependence on the choice of the tetrad We chose in section [3.2.21 to align the 
tetrad e° with the form (drj), which means that the corresponding observers follow 
worldlines always orthogonal to the constant time hypersurfaces. Alternatively, we 
could have chosen to align Gq with the vector (9,,), which would then in turn correspond 
to observers of constant spatial coordinates [9]. These two choices are related by a 
boost parameterized by the shift vector w^, and in the case of the Newtonian gauge, the 
parameter velocity which characterizes this boost reduces to This change in the 
tetrad used to decompose the momentum induces a transformation of the distribution 
function. Up to second order the two gauge invariant distribution functions (with 
obvious notation) are related by 

= Adi) _ (4-35) 

j{2) _p) , o^(2) , 

It can be checked that the equation satisfied by /^g^j is the same as ^^^^-f 
[equation (|4.2ip ] with dj^'^J^n'^n^ replaced by — $j^'n-', and thus our result is 
consistent with [15l[48] where I\^q ^ is used. 

What would be the best choice then? In our case (the same choice is performed in [3] ) , 
Cj = 0. We can check that if there is no scalar perturbation, then the coordinates of 
the acceleration defined by ai, = (e°)^Vpe° read 

a/ = (e'})'+He'}. (4.36) 

We thus conclude that the choice that we made corresponds to picking observers 
which accelerate the less, that is which accelerate only because of gradients in the 
gravitational potential, and thus this class of observers is the closest to freely falling 
observers. 



5. The gauge invariant Liouville equation for baryons 

5.1. The background and first order Liouville operator 

We follow the same method as for radiation. The main difference lies in the fact that 
n' is not a unit vector anymore. Consequently we cannot drop terms like UiU^ = f? . 
Furthermore it is not interesting to split the dependence in the momentum into a 
dependence in its energy and in the direction of its spatial part, since the magnitude of 
the spatial momentum (A) cannot be identified with the energy q° . At the background 
level, the evolution of 5°, (f and A are dictated by 



and we thus obtain for the background Liouville operator 

In order to compute the Liouville operator at first order, we need the evolution of 
the particle spatial momentum magnitude A at first order, which from equation (j4.3p 
is given in the Newtonian gauge by 

(1) 1 /A.o\{l) 



dxy'> _ 1 /dg°Y 

drj J 13 \ drj J 



A 



(5.3) 
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and this leads to the foUowing first order gauge invariant LiouviUe operator 

1 



/?2 



5.2. The second order LiouviUe operator 



(dxy^ _ 1 /dg; 



X 



At second order, the evolution of A is also obtained from equation (|4.3p . and in the 
Newtonian gauge we obtain 

(2) -I /^^o\ (2) i 

+2{^ -^)n'di^ + . (5.5) 

We also need to use the first order evolution of which in the Newtonian gauge reads 
do*\''' 

— ^ci^' -q°d^'^-q°{l3^-i'^ -n'n^)dj^. (5.6) 
ar] J 

Finally we need the first order photon trajectory given in the Newtonian gauge by 



dx^ 



i\ (1) 



(1) 



(5.7) 



With this we obtain the second order gauge invariant LiouviUe operator for massive 
particles 



drj 



(5.8) 



+ 



(2) _ H-(2)' 
IJ 



n n-' 



dq' ^ 



(5.9) 



.|5 2^ .^,_2|,i#(i)j^(i) 
dq'- dq' ^ 



5.3. The fluid limit 

If we choose a tetrad, not necessarily adapted to this bulk velocity, then the 
components of the stress-energy tensor are given in function of the energy density 
and pressure by 

^3„k 

(5.10) 



T°°^p + {p + P) [{u^f -1]=J : 
T°* = pu°u' 



(2^) 



T'l = pu'u' + PS'^ = 



(5.11) 
(5.12) 



In general for a distribution function g{q^), its associated stress-energy tensor 
conservation tensor V^T'"' is obtained by [T] 

aV,T^' = I L*[9]q'-^, (5.13) 
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where we recall that a is the scale factor. We can then perform a perturbative 
expansion of this expression according to the expansion (|2.7|) and obtain for baryons 
the perturbation of the continuity equation 



d3 



^ q°L*[g] ^ {T°°)' + mT°° + HTl 



(27r)3 



i2n) 

rjiOo{l) 



■oo{l) 



nr. 



(5.14) 
(5.15) 



a(v^r^°)^^^ 



(27r)3 



(5.16) 



jioo(2) j _|_ 3-^7100(2) _^ -^y'^(2) _j_ Q^j,oi(2) _ x|,'(2) (37^00 _j_ jii^ 

2(^ + $)a/f - 2*' (^3f ""(i) + f 'i^^') + 4r°*(i)a/ (4 - 

4 (31"°° + ^'^ 



2$ 



■00(1) 



3-^yoo(i) ^ _^ a/T'°*(i) - ^' (3T°° + T^) 



and the perturbation of the Euler equation 



(1) 



(27r)3 



(5.17) 
(5.18) 



a (V^T^' 



,(2) 



d3g* 



.2 .3-'?'^^<^^[^,^^^^5(^^] 



2$ 



4HT' 



oj(l) 



(5.19) 



(2^ - 4$)a,7$(i) + 



+ 2 (a^#f J;^^' - 39./#r^(^) 

In order to further specify the continuity and Euler equations, we need to 
express the perturbations of the stress-energy tensor components in function of the 
perturbations of the fluid quantities. Due to the strong coupling between electrons 
and protons through Coulomb interactions when matter is ionized, the baryons cannot 
develop significant anisotropic stress before recombination and thus the corresponding 
stress-energy tensor is of the form 

TI"' ^{P + p)u''u'' + Pg'"'. (5.20) 
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Indeed, the strong interactions between electrons and protons ensure that they follow 
a Fermi-Dirac distribution of energies in their rest frame (that is where they have no 
bulk velocity). We also note that since = 0, then the normalization condition of 
MqU" = —1 implies 

u°(i)=0, u"'-^^ =u'^^^u,(^iy (5.21) 

Furthermore we can identify the velocity and the spatial momentum up to second 
order in perturbation since 

u' = ^ ^ v' . (5.22) 

We then deduce the following perturbative expansion of the stress-energy tensor in 
function of the fluid quantities 

T°° =p, (5.23) 
^00(1) ^ p(i) _ (5.24) 

T°°(2) = pi2) ^ 2{p + P)u'(i)u'/^ , (5.25) 



T°' =0, (5.26) 

7^0^(1) = + , (5.27) 

J,o^i2) = (^p + p) + 2(/i) + F(i))u'(i) , (5.28) 

f^J ^ pj^3 ^ (5.29) 

^ ^ (5.30) 

T'^^^^ = _^ 2(p + P)u*(i)it^'(i) . (5.31) 

By using these expressions, we obtain the following gauge invariant conservation and 
Euler equations from equations (|5. 14115. 19)) 



a(V^r'^°) =p +37^(p + P), (5.32) 

a(V^Tn^'^ =p(i)'+3w(p(i)+p(i)) 

+ (p + P)ajM*(i)-3(p + P)*', (5.33) 

+ {p + P)die^^^ - 3{p + P)*(2)' + 2{p' + P')u'^^^u[^^ 

+ 2n{p + p)e^'^u[^^ + 2di [(p(i) + p(i))u'(i)] - 6*' + p(i)) 

+ 2($ + *)(p + P)diu'^^'> + 4(p + P)it*(i)9/($ - *) 

- 12**'(p + P)+ A{p + P)u'^^^u[^^' ~ 2$a (V^T^°)'^^ , (5.34) 
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(5.35) 



>(2) 



AH 



m 



2(4 + 'I') <^ (P + P) 



AH 



(p + P) 



4 (p + P) $9^$ + 29^$ (^p'l' + P(i' j - 8i/'{p + P)w*(i) 
2*a(VpT^^)*^^ . 



(5.36) 



5.4- Comparison with the coordinates approach 

In the literature, the fluid equation are usuaUy derived starting directly from the stress- 
energy tensor and in components associated with coordinates. It is straightforward to 
compare with our results since 



At first order we have 
and at second order 

a(v^T'"'^)'^' = [(v^r^°)(^^ + 2$(i) (v^r^°)^^^ 



(2) 



v(2) 



(5.37) 

(5.38) 

(5.39) 
(5.40) 



The presentation that we have chosen in the previous section makes the use of these 
relations simple. However, in order to succeed this comparison it should be noticed 
that HSl 



„H2) ^ y/(2) _ 2*1^^(1) , 



(5.41) 
(5.42) 

since the results in the literature for the fluid approximation are so far expressed in 
function oi . 
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6. Collision term 



6.1. General expression 

The general expression of the collision term is given by [701 [ZHl IHOI 
Cu.Jp) - ^u:^u^j27r)4 JdciJ Dp' J Bq'Sj, (p' + q - p - q') 
{ 9c{ci)UAp')M''Up',ci;p,q') 
- .gc(q') [V''^U^P{P)UAP)] M^'-^Ap, q'; p', q) 
+ [,9c(q) - .9c(q')] UAp')M^"'.riP' .<i;P,^')^^''\p{p)} , 



(6.1) 



where we used the notation Dp = S]j{p.p — m?) ^^^yi . This expression accounts for 
the incoming transitions in the second line and the outgoing transitions in the third 
line. We have neglected the Pauli blocking terms {1 — gc/2) since for the physics of the 
CMB, the electrons are much more diluted than the photons. However we include the 
stimulated emission in the fourth line and discuss later its relevance, and this requires 
the use of the stimulated emission matrix 



(6.2) 



In this definition we recall that the projectors S'q"^ are taken with respect to an 



observer having a velocity q/m, which means that SA'^qp = 0, and f"'' is the 



complex conjugate of f'^'^ . The stimulated emission in the case of polarized light 
is discussed in [531 [HHl ISI]- These treatments give the stimulated emission matrix 
in terms of the Stokes parameters. However, we do not want to refer to them 
since they imply and choice of the angle in the polarization plane basis, but it 
can be cheeked that the expression given here for the stimulated emission tensor is 
equivalent to the expression of these references when expressed in terms of the Stokes 
parameters. These parameters are an intermediary step between the polarization 
tensor and the multipoles which is not compulsory, and we thus intentionally bypass 
them in our derivations. In the end, this method will prove more powerful because it 
is computationally more straightforward. 

The Dirac function in equation (j6.ip ensures momentum conservation and 
A4'^'^^^{p, q; p', q') is the transition tensor for the process 

7(P) + e"(q) < > 7(P') + e"(q') ■ 

Its expression is given by 



A^^''„.(p,q;p',q') 



(6.3) 
(6.4) 



where 



p q ^ P:q 2 

p.q p'.q 



O 



Q^':0 ^ s^^'.s^^s'^lX^ ^ s^^.^{p)s^^Ap)s^Ap')s^%{p') ■ (6.5) 

The detail of its derivation can be found in [79j . Intuitively, it consists in transforming 
the incoming polarization tensor to the rest frame of the incoming electron, and 
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applying the Klcin-Nishina formula for the amplitude of the scattering process. 
Thanks to the reversibility of the scattering process, the transition tensor satisfies 

>('^%(p',q';p,q) = ^q'^.^q'';A^"".A(p,q;p',q')5q';5q/^, (6.6) 

as it can be checked directly on its expression (|6.4p . It is then possible to recast the 
collision term as 

CuUp) = ^u^^u^ l^qj Dp' J Dq' {2TTf6t, (p' + q - p - q') (6.7) 

X ge(q)^.(p')-^'^';.(p',q;p,q') {s^JJ, + M^^^^ip)} 

- ^u^^u^ J DqDp'Dq' {2tt)^S% (p + q - p' - q') 

X.ge(q){?7"^C/a/j(p)/M.(p)+AAq^,„^(p)r^(p')}x^';.(p,q;p',q'). 

The remaining two screen projectors at the very beginning of equations (|6.1l6.7p en- 
sure that the result is read in the required frame (usually the eosmological frame). 
The form of the expressions (|6.1l6.7p is manifestly covariant. First, the measures Dp 
and Dq, the Dirac function and the electron distribution function .gc(p) ^-rc Lorcntz 
invariant. Second, given the transformation rule p.46p for the polarization tensor, the 
screen projectors in the transition tensor and the stimulated emission tensor ensure 
that UAp')M'"'^0ip\ q; p, q') and /^,(p')X^';,(p', q; p, q')/r"(P) ^^'^ independent 
from the frame used to evaluate /^i^. Finally we note that the screen projectors 
'S'uk'S'u^ are defined with respect to the observer having a velocity u, which is the 
velocity with respect to which the collision tensor in equations (|6.1l6.7p is defined (in 
our case = (dyy)^). Hence they ensure that the expression of the collision tensor 
is independent from the observer (that is the frame) used to define the unit polariza- 
tion tensor Uap appearing in the third line of equation (|6.ip and the fourth line of 
equation (j6.7p . and also independent from the observer used to define the polarization 
tensor fap appearing in the stimulated emission tensor. Additionally these screen 
projectors makes it straightforward to check that the collision tensor transforms as in 
the rule (fOg]) . 

Given this discussion, it is in principle possible to express this collision term in 
the fundamental tetrad basis, in order to specify further the collision tensor. How- 
ever, all the screen projectors appearing in the transition matrix and in the stimulated 
emission tensor are taken with respect to the observer with velocity q which varies 
throughout the integration J Dq. In order to simplify the computations, it would thus 
be helpful to use an intermediate tetrad basis adapted to each momentum q in this 
integration. Additionally, instead of expressing the collision tensor in the fundamental 
frame or tetrad, one would prefer to compute it, that is to obtain it explicitly in func- 
tion of the radiation multipoles, in the electrons rest frame. In a second step, we would 
then transform the collision term in the fundamental frame using its transformation 
rule which is given by equation (|1.46p . and the transformation rules of the radiation 
multipoles and momentum components under this change of frame that we gave in 
section 11.71 We thus need for this final step to determine up to which order in the 
electrons velocity we need to compute this change of frame to be consistent with the 
perturbative expansion of the Boltzmann equation. 
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6.2. A double perturbative expansion 

When performing a perturbative expansion, we need to expand both in the metric 
perturbations and in the vefocity perturbations. However these perturbations are 
of different type and also have a different magnitude. The magnitude of the met- 
ric perturbations are typicaUy of order 10~^ for the processes involved in the CMB, 
and thus the bulk velocities of the radiation and matter distribution are also of the 
same amplitude. However the temperature of radiation (Tr) and the temperature of 
electrons (To), which are nearly identical around recombination arc approximately of 
order Ti/too — Tc/ruc ~ 10~^ and thus the thermal velocity of electrons is of order 
10~^. Consequently, when going up to second order in metric perturbations, we would 
need in principle to keep perturbations up to fourth order in velocities as they can 
in the end lead to terms of order {Tj./meY, T^Tg/ml and (Te/mc)^, which can be 
comparable to the second order perturbations of the metric if the numeric coefficients 
in front of these terms are quite large. It can be shown [521 [521 [HI] that in fact, due 
to Comptonization processes, these terms will be of order (T/me)^ x (Tj — Tc)/T and 
will be thus completely negligible in our context. However we would certainly need to 
retain terms up to third order in the baryons velocity, that we will call accordingly to 
this discussion of order 3/2 in the metric perturbation (though there is no such half 
integer perturbation for the metric). 

We will now apply this method to perform a perturbative expansion of the collision 
term 

• in section 16.31 for electrons having no thermal dispersion and choosing for this 
computation the rest frame of electrons which is also the rest frame of baryons, 

• in section [6.41 for electrons having a thermal dispersion, choosing the same frame 
for the computation. 

• in section 16.51 for the general case, accounting for the bulk velocity of electrons 
by transforming the result obtained in the baryons rest frame toward the 
fundamental frame. 

6. 3. Distribution of cold electrons with no bulk velocity 

It will prove convenient to calculate the collision term in the Boltzmann equation 
starting from a very simple case. We will first assume that the free electrons have 
only a bulk velocity and no thermal dispersion in their velocity distribution function 
go{q^). Furthermore, we will choose to align the first vector of the tetrad Gq (an 
adapted tetrad) with this bulk velocity that is to work in the baryons rest frame in 
order to specify the collision tensor given in equation (j6.7p . In that case the electron 
distribution function is explicitly given by gc{<t) ~ (27r)'^(5£)(g*)no which makes the 
integration on electrons momenta trivial. In the baryons rest frame, the velocity of the 
incoming electron is aligned with the tetrad Gq, and we obtain from the momentum 
conservation in the collision 



The second term in the brackets of this expression is of order (Tr/me)^, that 
is comparable to second order perturbations in the metric, and according to the 
discussion in section 16.21 they can be ignored. We then need to expand the Dirac 



1 



O 




(6.8) 
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function. Wc first integrate on the spatial components of the outgoing electron 
momentum since the integrand does not depend on it in the form (j6.7p of the collision 
tensor. The remaining one-dimensional Dirac function is handled by following the 
standard steps that can be found in [5], and we obtain 



2rrir, 



dp' 



-Mp" 



P 



^Dip'" + 9'° - 9°) = 5],[p'^-p'") + ^2!!!' ^'"^ ^o^'^^P'^-P") ■ (6-10) 
Additionally, using the fact the tetrad corresponds to the rest frame of all electrons, 

(6.11) 



— = — + (1 - n'n[) , 

p° p'° Too ^ ' 



which is more conveniently expressed as 



(1 - nVQ 



(6.12) 



P° mc 

As a consequence, the term in square brackets in the definition (|6.4p of the 
transition matrix is of order {p° /nie)^ ~ {Tr/rrieY and is comparable to second order 
perturbations of the metric. This would not have been the case if the tetrad was not 
adapted to the incoming electron velocity and in that case the term in brackets in 
equation (|6.4p would be of order Tr/me- This property is important since, given the 
discussion in section [6.21 they can be ignored. Finally using 



the collision term is expressed as 

Cab {p") = n,aT [Cab,T (p'*) + Cabjl (p'*)] 

where 



Cah,T(/) =P° 



in 



S^Xfcd (p") - fab [p'') 



CabMp'') = I dp'V 



2o' Vf^o\2 



d^' 

47r 



/„.(p'^)-/(/)^ 



(p° ) 2 + (y°)2 - 2p'yn'n[ 
2mo 

[1 + lip'^)] SlStfcd {P 



d6]j{p"'-p°) 



(6.13) 
(6.14) 
(6.15) 
(6.16) 



Qp/o 

^2^P-^fab{p^). 



We recall that in these expressions, the screen projector 5'^'' is defined with respect 
to the photon with momentum p, and thus S'^'^ f cd{p'^) 7^ I{p'^)- Note also that 
S'^'^S'^^ = 1 + (n.n')^ [79], and this helps recovering from the trace of the Thomson 
term (j6.15p the standard form of the Thomson collision term for unpolarized radiation. 
Using the integrals p.52p we find explicitly the Thomson and the recoil terms which 
are given by 



Cab,T{p'') „ 1^ 



-/(p°, n') + J(,(p°) + l/,<,(p°)n'=n'^ - ^i?ed(p'')n'=n'^ 
10 5 



-Pabip^^n') - ^Iab{p°) + \Eab{p°) 

10 5 



TT 



+ -j^'^^abc'n 



^V{p°y) + \va{p")n'' 



(6.17) 
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hip") - + ^Icdip'^Kn'^ - ^Icdkip^Kn'^n'^ 

5 10 70 

'^E,d{p")n^n'' + l-E,dH{p°)n'n^n'' 



|p,fc(p°,n^) [2 + 4ip") + 2/0(p°)] + [1 + kip")] (2 



^Eabip") - lEabc{p''K - iBUpl^b^na 
5 7 5 



TT 



iie„fc,n^|y(p°,nO [2 + /^b") + 2/0(p°)] + [1 + (s+p"^) 

-\V^{P") + - II . (6.18) 

We recall that, according to the notation of section fT.7.21 I'af {p°) = P°-^Iae{.P°)- We 
see clearly by examining the antisymmetric part of the Thomson and recoil terms, i.e. 
the terms proportional to icabc'Ti'^i that the circular polarization is not excited and thus 
remains null if it is initially so. From now on, we will discard the antisymmetric term 
in the polarization tensor and in the Boltzmann equation. Additionally, since there is 
a small amount of energy transferred from the electrons to the photons in the recoil 
term Cah, r, there is also energy transferred from the photons to the electrons, and 
it is incompatible with a cold distribution of electrons. In order to have a consistent 
treatment of the collision term, we can still work in the electrons rest frame but we 
have at least to consider a thermal distribution of electrons. Consequently we turn 
to this case in the next section, but the results obtained for a cold distribution of 
electrons can be used to derive the collision tensor for this more general case. 



6.4^. Thermal distribution of electrons with no bulk velocity 

If the distribution of the velocities of free electrons is thermal then we can apply the 
previous method for each of the electrons. That is for each momentum q in the integral 
on distribution function .gc(q), we choose a tetrad field Gq whose timelike vector Go is 
aligned with this momentum (gq = q/rrio). We then calculate the collision rate per 
electron defined in equations (|6.15p and (|6.16p in this frame [that is (Tt^cC*-^ t(^''^) 

and a^ficC^i r(p'^)]i f-nd since these quantities transform similarly to equation p.46p . 
then the collision term is expressed by 



Cabip" 



P 



C 



(6.19) 



(27r)3g° 

In practice this means that we need to use the transformation rules of the radiation 
multipoles, of the energy and of the momentum direction, in order to express 
everything inside this integral in function of the multipoles taken in the fundamental 
frame. The interest lies in the fact that integrals on the electron distribution function 
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with odd powers of v (the velocity of the electron considered) vanish. Since we need 
only to go up to third order in the thermal velocities, then in that case we can restrict 
to second order in thermal velocity. 

However we have many reasons to follow a much simpler track. First the terms 
coming from the second order change of frame (that is second order in v) will 
lead after integration on the electron momentum to terms having a factor at least 
Te/nie = (T + dTe)/me, and following the discussion of section [621 this factor can be 
comparable to first order perturbations in the metric {T /me) or second order {STe/rric). 
Since we stop our expansion at second order, we can have in each of this term either 
a quantity whose smallest non- vanishing perturbation is the first order (lag for £ > 1) 
or a quantity whose smallest non- vanishing perturbation is the background (/o). If 
it is the latter case, then the term contributes to the famous Kompaneets collision 
term |85) . If it is the former case, then the term contributes to what can considered as 
corrections to the Kompaneets term for anisotropic radiation, but since it is linear in 
first order perturbations, this type of term cannot contribute to the bispectrum |44| . 
Similarly, the terms coming from the recoil term have an overall factor of p° /rUe 
which is typically of order T^/mc, and for these terms the same reasoning applies. 
Consequently, when interested in the bispectrum, we should in principle just keep the 
Kompaneets term. However the expression of this famous Kompaneets collision term 
which is contained in the total collision term (|6.19p is 



(6.20) 



and it can be checked that for I(d{p°)/2 following a Bose-Einstein distribution of 
temperature ~ To, then the Kompaneets collision term exactly vanishes, and if 
these two temperatures are close, as it is the case when recombination occurs, then it 
is of order {Tr~Tc)/mc. The Kompaneets contribution, which is at least a second order 
quantity in the sense of the discussion in section [621 bears thus only a linear response 
in the (primordial) first order metric perturbations through T/me(r/^-' — Tc^^)/T. 
Consequently it cannot contribute to the bispectrum generated by evolution, and we 
will ignore it in the rest of this paper. Note also that on the form ()6.20p of the 
Kompaneets collision term, it is obvious that ^ = aCcd, k /p° does conserve the 
photons number density [2l [46] and not the energy density (it is known to induce 
spectral distortions) as it is not the case on the expression (4.38) given in [48j. Hence, 
the integrated counterpart of the Kompaneets term, C^, does not vanish and it is only 
because we focus on the bispectrum that we may discard it. 

Gathering all these remarks, we conclude that for a thermal distribution of electrons 
with no bulk velocity, the collision term is given by 

Cab{p') = n,cTTS^,S^CrAMpl + 0(1) [0{T,/m,) + 0(p7me)] , (6.21) 

and that the terms of order O (To/me) and 0{p°/me) are irrelevant for the purpose 
of computing the bispectrum generated by evolution. This method which consists in 
working in the baryons rest frame simplifies the physical interpretation of the collision 
tensor. Indeed on the form (|6.17p and (|6.2ip we clearly see that in this frame in the 
tight coupling limit (that is when Cab. t 0) there is no polarization {Pab = 0) and 
the radiation is moving with the baryons as if it was a single fluid since the only non- 
vanishing multipole is X0. This terminates rigorously the discussion in section II. C. 2 
of [Hj (see also the discussion at the end of section [^^ . 
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If the distribution of electrons has indeed a bulk velocity in the cosmological frame, 
then we just need to transform the result obtained in the baryons rest frame to the 
cosmological frame, and in this process, we keep only the terms which can contribute 
up to second order when performing perturbations. For the sake of comparison with 
literature, we only report the intensity part of the collision term and present here 
rather than C. We obtain 



(1 - n.v) 



(6.22) 



- l/^,(p°)n^n^n.v - l/e(p°)n=n.v + l/^(p°)n^n.v 

- /;(p°)v.v + ^/;'(p°)v.v + ^I'^ip") (n.v)^ + /;(p°) (n.v)^ 

- ^E,a{p'')n'7i''{l + v.n) + ^EUpln^n^^.n + ^E,d{p°)n'v^ 

-\BUP°y^'"'vhnfn- + ^B'^^{p°)tf''<'vhnfn' 



0{Tc/mc,p°/me) , 



where r' = anoCr-il^l. 

Though we have also obtained the full collision term for the polarization part, wc only 
report here the result of its energy integrated counterpart, in order to simplify the 
expressions obtained. 

1 



T-Sab 



-X(n') +T0 + Ij.^n^n'^^ (1 - n.v) 



(6.23) 



1, 

2' 
1, 

5 
3 
5 



^ 'T r , '-7- r fl 

-Xcn n.v H — Xcrfn ?? n.v 
2 5 



Icdn'^v'^ + 42'0n.v - Xgv.v + 7^0 (v.n) 



-Ecd'n'n'^il + 5v.n) + -Ecdn'^v'^ - -Bcd^' '■"VhUfn" 



1 3 

Y^^ab + -^Sab ] (1 +3v.n) 



TT 



+ ^£acn''vb - ^lacn^Vb + ^B^^'^ebhdv'^ + ^laVb - T(t)VaVb 

+ 0{Tjm,,p"/m,). 

§ Note here that rio is the electron number density as seen in the electrons rest frame. Though 
we have changed the frame for expressing the collision term we keep this definition for the electron 
number density since this is the most physical. However, a change of frame only affects the number 
density by a factor 7 ~ l + v.v/2, and since the background collision term vanishes, the second order 
perturbation of the collision term is not affected by the precise choice of frame used to define the 
electron number density. 
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The trace part of this energy-integrated coUision term as weh as its non-energy 
integrated counterpart which is given in equation (|6.22p reproduce the results obtained 
previously in [ISl H51 [55] for unpolarized radiation. 



6.6. Multipole decomposition 

This last expression can then be decomposed in multipoles, in order to separate the 
intensity from the linear polarization part. 



je + i) 

(21 + 3) 



(6.24) 



5} [3X0 Wai] + (5| 
1 „ 



10 



(^-l)(£ + 3)^ , 



+ 6j 



^-p-|-efce(o,f''ga,_^) 



(6.25) 



10 



-C*[BU = - B 

T — 



(^-l)(£ + 3) 



^ (2^4-3)(^-hl) 



-^BaebV + g(at_it'of) 



(6.26) 



^+1 



^bc(a 



4 6cc ^ 6-7-r 



We recall that the multipoles are projected so all indices a,b,c. . . in the above expres- 
sions could in fact be replaced by z, j, fc . . .. 

In the case where there is no polarization, the expression (|6.24p for the intensity part 
of the collision is nearly consistent with equations (64-67) of [57] . Indeed the only dis- 
crepancy is the coefficient in front of I^j^vi^ai'^a-i) i which in our case is 7 and in the case 
of [57] is 3. In their notation, this discrepancy arises from a missing factor 6 in their 
equation (60) in front of the factor {v%ecf' . The terms in the first bracket of the r.h.s 
in this equation can be traced directly to the expansion of [7b(1 — v%ec)]~^ , arising 
from the change of frame between the baryons rest frame and the fundamental frame. 
Indeed, in the definition 14.61 apphed to the colhsion tensor, C[] is a scalar but (noting 

here the tetrad in the baryons rest frame) [p°)^(lp° ~ {p°)^dp° [7(1 — ri^Vi)\ , and 
in the definition 14.11 of there is an additional factor \/p° = 7(1 — n^Vij/p". As 
a consequence, the factor —3 in front of the term pjiVgv''g in equation (63) of [57] 
should be —7. 



In the tight coupled limit, that is when we can neglect (when compared to the 
Liouville operator) the r.h.s of equations (|6. 24116. 26]) . then it is obvious that 

= B,, = X,, = 0, if £ > 3 , (6.27) 



p. . — R. . — n 
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I I = 4I([iVi . 

We recover of course that since there is no polarization in the baryons rest frame, as 
noted at the end of section 16.41 then there is no polarization in any frame, and we 
check from the comparison with a perturbed stress-energy tensor with no anisotropic 
stress given in equations (|5. 28115. 3T|) that the multipoles Xa and Tab, which vanish in 
the baryons rest frame, are non-vanishing only because of the change of frame, and 
consequently in this tight coupled limit there is a quadrupole but no anisotropic stress 



6. 7. Perturbative expansion of the collision term 



Though this result might appear to be general since we have not performed a 
perturbative expansion around a FL space-time, it is not general since we have kept 
terms quadratic in the change of frame velocity v only when they multiply Xg, since 
this is the only quantity which will have a background contribution. The full second 
order transformation (in v) can be computed using the transformation rules at second 
order in v for the multipoles involved in the Thomson collision term (|6.17p and we have 
reported them in [Appendix A[ This full result at second order in the bulk velocity 
V has already been derived in [SHj , but without polarization. By developing with the 
expansion (|2.4p the intensity, the electric and the magnetic multipoles along with r' 
and V, we can perform a perturbative expansion of the collision term following the 
same notations as for the Liouville operator expansion, that is as in equation (j4.5p . The 
background collision term vanishes. Since we have decided to discard the first order 
magnetic multipole of radiation (because we can neglect first order vector and tensor 
perturbations), there is no contribution to the magnetic multipole of the collision term. 
The first order intensity and electric multipoles read 
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At second order, the linear contribution is given by 
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As for the quadratic contribution, it is given by 
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6. 8. The fluid equations for baryons 

The total stress-energy tensor is conserved 

V^'^bTr = ^M^b" + = . (6.36) 

This arises from the Bianchi identities, and thus the action of baryons on photons 
is opposite to the action of photons on baryons. This means that we can define the 
force resulting from the action of photons on baryons and the force resulting from the 
action of baryons on photons by 

V^T^^"- = F^^^ = -V^Tr = -K^r • (6.37) 
The expression of the force in the tetrad basis is further given by 

and from equations (j5.13|) and ()1.59p , its components are easily related to the moments 
of the radiation collision term by 

aF;^,= -C*[I]\ (6.39) 

aF:^^= -ic#[/r. (6.40) 

It is then straightforward to deduce the perturbative expansion of the force, which 
is inherited directly from the perturbative expansion of the radiation collision term. 
Note that in particular F°}^^ = 0. 



7. From PSTF multipoles to normal modes components 
7.1. The normal modes in Fourier space 

The equations (|4.25II4.2 j)) . (|4. 31114. 3^ . (|5. 34115. 36|) with equations (|6. 30116. 35)) and 
(|6.39p are a key result of this paper, since they provide the full second order Boltzmann 
hierarchy for the coupled system of baryons and photons. However in order to 
make contact with the standard approach of CMB numerical integration, it is more 
convenient to describe the angular dependence of the radiation functions using the 
normal modes components |15] . And in order to perform a numerical integration in 
the most simple way, we also turn to Fourier space in order to solve only for a time 
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integral. It has been shown [55] that we can convert the PSTF multipoles into normal 
modes components by decomposing the quantities X, £ and B according to 

X{x^y) = J2 / —-^Xr{k,v)Gf^{k,x',rin = ^ , (7.1) 

where X stands for X, £ or B. Here we have used the notation k, that is a boldface 
font, for a mode in the Fourier transform though it is only a three dimensional vector. 
The modes G^(n*) according to which we decompose are defined by 

GL(k,x^n) =xi-e''=^^V^™(n), (7.2) 
GL(k, n) = -GL(k, x^ n) = I-^e''^^- Y^"(n) , (7.3) 

with 



The components (k, 77) are the time-dependent components in Fourier space of the 
spherical harmonics decomposition, also called normal modes. Note that we have also 
factorized the background intensity so these quantities are now dimensionless, but we 
did not divide by 4 as is usually done when we want to have a quantity which can 
be interpreted as a temperature. Instead we decompose I /I which is the fractional 
variation of the energy density of radiation as it is what is really measured. The 
A£m(k, 77) can be obtained from the PSTF multipoles using the formulas 

Xr(k,r/) = iiV.A, j ^^e-'^''^'y;^iajx-^) , (7.5) 



The functions y^-^ are null if one of the indices in at is o, since it is already the case for 
the indices in the intensity, the electric and the magnetic multipoles. These quantities 
are built in details in |i69, , and in the case where to > they are defined by 



yfm = (-ir 



21+1 (£- m)! 



An (^ + to)! 



1/2 [('?-™)/2] 



X ((5f ^ + i4'') . . . + i^") S'f^' . . . sl'-^' 5''-'^+'''-^^+^ . . . 5''-''' 

If TO < 0, then it is defined by = (-l)™3^^'^jf„. The PSTF muhipoles can be 
conversely obtained from the normal modes components by the inverse formulae 



-^(e^-ie^). (7.13) 
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- -i^, I ty'^mK ,) - (7.11) 

7.2. The basis for the decomposition of tensors 

In order to decompose the metric perturbations and the fluid perturbations according 
to this spherical harmonics scheme, we start from the vectors 

e(o) = -e3, 6(1) = -i= (ei +ie2), e(_i) = --i= (§1 - iea) . (7.12) 

These are not true vectors since they hve on the background space-time and their 
indices are raised and lowered by the euclidian metric Su and its inverse 5^'^ . Their 
associated forms are thus 

e(0) = -e3, e(i) = ^(ei+ie2), e^-^) = -1^ 
^/2 ^ ' V2 

We then build a tensor basis in Fourier space out of this basis. We start from the 
scalar basis which in Fourier space is given by 

Q("' =exp(e'^'^^^') , (7.14) 

and we align k with §3. It is then used to build the scalar basis of higher rank tensors 
according to 

Qf^ -9,QW/fc = ie-(?)Q(°), (7.15) 
Qf] ^ a(,9,)Q(°)/fc2 ^ (_g(o)e-(o) + 5jj/i) , (7.16) 

since perturbation variables "live" on the background space-time. We then remark 
that 

Q(")=Gjo, n'QW=Gfo, nVQ^°) = ^Gfp , (7.17) 

and this enables us to make contact with the normal modes decomposition first 
performed for the directional dependence of radiation. For vectors, we follow the 
same method and we start from the basis 

Q(±^)^ie-(f)QW, (7.18) 

from which we can build a basis for vector type perturbations of higher order rank 
tensors (we restrict to rank two tensors, since in the problem at hand we have no rank 
higher than 2) 

Q't,''^-d,iQf'lk = i^^,QfK (7.19) 
We then make contact with the normal modes decomposition by noting that 

= Gf±i, n^n^Q^^j"^ = 4=Gf±i . (7.20) 



Note the difference in the definition of Q^^ in equation (|7.18p with respect to the 
expression given in [15j . However it agrees with the expression given in [90j up to 
factors (—1)'" which arise because of a different convention in the spherical harmonics. 
Finally for tensors we use the basis 

Q'u' - -\/!e-re-f , (7.21) 
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which conveniently satisfies 

nVJQ±2 ^ gi^^ (7 22) 

The (second order) vector perturbation of the metric is then decomposed as 

^f\x-') = J E Qi^'^Hml^JiKv)- (7.23) 

^ m— ±1 

Similarly the (second order) tensor perturbations are decomposed according to 

= / jf^ E Q?.r\^)H'r^HKri) ■ (7.24) 

As for gradients of the first order and second order scalar perturbations ($ and 5") 
they are decomposed similarly to the vector perturbations in 

dM^^) = / JYYT^ E Q/^'"^(k)fc('"^<I>(k,ry). (7.25) 

^ ^' m=-l 

This is equivalent to decomposing the Fourier mode k into a scalar part (in the 
direction ejg)) and a vector part (in the plane orthogonal to e(o)) according to 

■m— — 1 

The relation of these components to the Cartesian components of k are 

fc("' = -fc3, fed) -L(fci-ifc2), fc(-i) = --L(fci+ifc2).(7.27) 

For the second order perturbations (or for first order perturbations when we solve the 
first order equations) we choose to align the mode k with §(3), which implies in the 
previous decomposition that fc*^^-* = fc'"-'^-' = 0, and fc'"-' = — fc. However for first order 
quantities appearing in quadratic terms in the second order equations, this is not possi- 
ble because of the convolution generated by the Fourier transform. In the following we 
will also use the notatior[ij] k = k/fc with the same decomposition as in equation (|7.26p . 

As for the second order velocity of the electrons v, it has a scalar degree of freedom 
and a vector degree of freedom according to the decomposition (|2.6p . and the Fourier 
components of this vector degree of freedom could be decomposed similarly to 
that is according to equation ()7.23p with m running from —1 to 1. However since the 
velocity of the electrons always appear with tetrad basis components, it is not very 
convenient to go in the coordinate basis in order to split into scalar and vector modes 
and then transform back to tetrad components. We use instead for Va the same type 
of decomposition as we did for Xa ■ Namely, this decomposition reads 

-^i-^) -jrj ^2-'""'' ily^-'-^^^^r^)- (7-28) 



II Note here that the hat denotes a unit vector and has nothing to do with a gauge invariant variable. 
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7.3. Comment on the SVT decomposition 

We have arbitrarily chosen to perform the muhipole decomposition on z)" and not 
V^. If we use for this multipole decomposition, then Vq matches the scalar 
component and V±i match the vector components in the usual SVT decomposition of 
equation ()2.6|) . We did not follow this approach for the baryons velocity, but we did 
it for the vector degrees of freedom in the metric and also for the tensor degrees 
of freedom H^'^ . We see from equation (j5.4ip that ^ V'^ and consequently, when 
performing the multipole decomposition using -O" (which is what we have done in this 
paper) the normal modes components components arc different (vm 7^ Kn)- It would 
thus be misleading to name vq a scalar component and wj-i vector components. For 
the sake of the discussion here we will nonetheless use the terms coordinate SVT and 
tetrad SVT to name these different decompositions. Essentially the difference comes 
from the fact that when we perform the coordinate SVT decomposition of a vector, the 
decomposition is made on a Euclidian space-time which is the background space-time. 
Indeed the indices on the mode are lowered and raised by a simple kj = Sijk^ . The 
coordinate SVT decomposition of tensors is not motivated by physical reason but just 
by the Fourier transformation. On the other hand, in the tetrad SVT decomposition 
we use the tetrad indices i, that is we work in the locally Euclidian frame of the 
physical space-time, and indices arc lowered and raised by Vi = 6ijV^ . The fictitious 
background Euclidian frame is different from the physical locally Euclidian frame. 
As it can be seen from equation (|5.4ip . a local volume expansion higher than the 
average (vf) leads to a difference in what is considered as a length, and the SVT 
components found in one decomposition are a rescalcd version of those obtained in 
the other decomposition. If we had chosen not to neglect the first order gravitational 
waves, then the r.h.s of the second of the equations (|5.41[) would be supplemented 
by a term 2Ef^^V'''^^\ Again, this would lead to a difference in the direction used 
to perform the decomposition and the scalar part of a vector in one approach would 
mix the scalar and the vector components of the other approach. Note that since 
we have neglected the first order vector and tensor perturbations of the metric, the 
decomposition of the second order vector and tensor perturbations of the metric in 
either approach leads to the same result. A more detailed discussion on the issue of 
SVT decomposition can be found in [HT] . 



7.4. Transforming the Boltzmann equation in normal modes components 

Since the second order equations arc quadratic, we need to know how to compose the 
. This is deduced from the composition rules of the spherical harmonics, and the 
key result is 



E 

(*l+f2-*3) ' 



' {2ii + 1)(2^2 + l) ^£3(mi+m2)^^30 

47r(2£3 + 1) ^imi<?2m2 Wi0^20 



where the Ci^™'^. ^ are the Clebsch-Gordan coefficients. A set of useful relations can 
be obtained by considering the cases where £1 — £2., £1 
are reported in [Appendix D 



(7.29) 

as can 

2 - 1 or £1 = £2 + 1- These 
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7. 5. The Boltzmann hierarchy in normal modes components at first order 



Throughout this section we will use the following definitions that will simplify the 
notation 



2TAT 



1 7>'m 



(£2-m2)(^2„s2) 



{l±m){t±m + l){P - s^) 



2P 



. =V 



(7.30) 
(7.31) 
(7.32) 

■a 2 

At first order we choose to align the mode considered k with the direction §3, 
and the dependence in k becomes only a dependence in its magnitude k. The set of 
equations obtained at first order is (dropping the obvious dependence of all quantities 
in r] in the rest of this paper) 



\m _ '"- ±1 \m _ 1 

Ap = — , Afi — m- 



1 /(£ + m=Fl + l)(^-TO±l) 



£#(i)[X]?(fc) =X;°(fc) + fc 



2£ + 3 



2i- 



c*'-^^[i]1{k)^i'/^{k) + k 



2£ + 3 



(7.34) 
(7.35) 



C*^^\lf,{k) = t' {-Xe°(fc) + 5°i°{k) + A5lvo{k) 

im - V64°(fc)' 



c*^''>[£]Tik) = ? l-srik) - s',^ [x°(fc) - Veslik) 



(7.36) 
(7.37) 



As already mentionned, the first order magnetic modes are not excited since the first 
order vector and tensor modes are negligible, so we did not report them in the above 
equations. Additionally, the intensity and electric multipoles are only excited for 
m = 0, that is the reason why we also only reported this case. 



7.6. The Boltzmann hierarchy in normal modes components at second order 

At second order, a Fourier mode k on a quadratic term will appear as a convolution 
on modes ki and k2 whose sum is k. It implies an integral of the form 

/C(ki,k2,k)^ I ^!|l|^5|,(ki+k2-k), (7.38) 

that we abbreviate in JC. We can only align k with the direction §3, but not ki or k2 
at the same time. However for a first order quantity the components X™ for a mode 
in a given direction ki can be obtained by rotating the components X'^ obtained when 
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we had decided to align this direction with e^. Namely this rotation leads for a mode 
k to 



4tt 



^r(k) = ^ ^^Y*'"\k)X', ike,) , 

and in particular for the first order velocity 

{}„(k) = -fc(")7)o(fce3). 
Note that the scalar product of two mode vectors is given by 



(7.39) 
(7.40) 



ki.k2 = fc(°)fcf - fc(^'fc(-^) - k[~'hi'^ = J2 -(7.41) 



7.6.1. The Liouville operator We finally obtain for the linear terms in the second 
order Liouville operator 
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(7.44) 



I iv^^ 

2-1^^ Rm(2), 



2£- 
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As for the quadratic terms, we can use the composition rules given in [Appendix D| to 
obtain 
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(7.45) 
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( 1 n tft'™ 



(7.46) 



^ n=-l 

1 n I T^-rn 

Tl=-1 

X 



(-«) 



(^ + 2)A:<-"M vl/(fci) 



^4* 



[(^ + 3)A;i")|.(fci) + (4"^ + - l)fci"' ) 
'(fci)£r(k2) - <l>(fci)£;™(k2)} , 



^ ^ n— — 1 

This last expression and its version in PSTF multipoles (|4.34p . can be traced directly 
to the lensing term in the expression (j4.30p and we recover the well known result 
that at second order, gravitational lensing generates magnetic type multipoles out of 
electric type multipoles. 



7.6.2. The collision tensor Using a similar method applied to the collision term leads 
to 

C#(^)[X]r(fc) = r' [-if'\k) + 5KST\k) + 45}4^{k) + 5|P'"(2)(fc)] , (7.48) 
where P™'^^'(fc) is not null only if — 2 < m < 2 and is defined in that case by 

(7.49) 



p™(2)(fc) ^ _L \lf^\k) - ^/64"^'^(fc) 

For the electric and magnetic type collision terms, their second order linear 
components read 



(7.50) 
(7.51) 
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The quadratic terms are then given by 



1 n fj^rn. 
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X™-"(k2) - ^/64™-"(k2)J i)„(ki)| 
2/CT'^'^(ki)C#(i)[X(i)]^(k2) , (7.52) 



1 nf-ry-m 
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WnCkl) 
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(7.53) 



— - ^ "Arz)„(ki)4"-"(k2) 

^ ^ n— — 1 



(7.54) 



^ "A^'f}„(ki 



^4'"""(k2) - ^x/6X2"-"(k2) 



7. 7. Continuity and Euler equations for baryons 



As discussed in section 11.51 we describe the baryons by cold matter that is by a fluid 
with no pressure. From equations (|5.33p (|5.34p and (|6.37p we obtain the continuity 
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equations 



(7.55) 



ri=-l 
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/5(fc2 



6*(fci)^'(fc2) 



2^(_1)» $(fc2)-#(ft2) 4"^«-n(kl) 



n=-l 
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1 - 4 

-Xf (k2)w-„(ki) - -■()„(ki)z}_„(k2) 



where R= pjX. 

From equations (|6.37p (|5.35p and (|5.36p we obtain the Eulcr equation for baryons 



lf'\k)+4vl^>ik) 



(7.56) 
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—-t'K, 
3R 
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(7.57) 



7.iS. The road to numerical integration 

The set of equations which need to be integrated is very intricate. First the scalar 
(to = 0) first order equations need to be integrated for each Fourier mode magnitude 
k (see iminilE^). This encompasses the first order Einstein equations to determine 
the first order scalar potentials, the baryons continuity equation (j7.55p and Euler 
equation (|7.56[) . together with the Boltzmann infinite hierarchy of coupled equations 
for the intensity and electric type multipolcs given by the first order Liouville operator 
multipoles (|7.34p and (|7.35p equated with the first order collision multipolcs (|7.36p and 
(|7.37p . The first order Boltzmann hierarchy couples the modes i with modes £±1. The 
results obtained depend only on the magnitude of a given mode since we chose to align 
the mode with the azimuthal direction of the spherical harmonics. The multipoles for 
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any direction of the Fourier mode is obtained by the rotation (j7.39p . Finally, we obtain 
from this first order numerical integration a first order transfer function T"^'^^'^ (k, 77) 
defined by 

ir«(k,^) =r™«(k,,7)xrWC(fc), (7.58) 

where C(fc) is the primordial curvature perturbation in comoving gauge. 

At second order the numerical integration has to be performed in the same 
manner. We have to integrate simultaneously the second order Einstein equations 
in order to obtain the metric perturbations and ijjj (the expressions 

of the second order Einstein tensor can be found for instance in [331 [Ml [551 f55])- 
the second order continuity equation for baryons (|7.56p and the second order Euler 
equation (|7.57p . together with the second order Boltzmann infinite hierarchy for 
intensity, electric type, and magnetic type multipoles which is obtained from the 
second order Liouville operator multipoles (|7.42H7.47)) equated with the second order 
collision multipoles (|7. 48117. 5¥|) . After numerical integration we will obtain a second 
order transfer function 7^"'''^''(ki, k2) defined by 

\if^\k,ri) = /Cr;"(^)(ki,k2,ry)Xr(ry)C(fci)C(fc2) ■ (7.59) 

However, we already see that the Boltzmann hierarchy at second order is much more 
complex than its first order counterpart. First it depends on the two modes ki and 
k2 (though their sum k is chosen to be aligned with the azimuthal direction of the 
spherical harmonics) inside the convolution in quadratic terms. Second, it couples a 
mode (^, m) with the modes ± 1, m ± 1). And finally, one problem which needs to 
be solved is the determination of t'^^K Indeed, photons can only scatter off free 
electrons and its correct expression is r' = arioXoCTT, where is the fraction of 
ionized electrons. Both n"c^ and x^c^ contribute to the first order perturbation since 
f'*-^-* = acTTx'i'^fic + aaTXcfi^c^ ■ Though fSe^ is easily obtained from the first order 
equations using fii^^ /fie = pi^^ /pa, the perturbation of Xc requires to determine and 
integrate the equations which rule the physics of recombination on the first order 
perturbed space-time. 

Conclusion 

In this paper we have extended our previous investigation [45] to the case of polarized 
radiation, also including the Compton scattering of photons off free electrons, in order 
to obtain a fully consistent second order treatment of the radiation transfer. We 
also studied the case of massive particles in the kinetic theory and checked that 
in the case T/m <C 1 it was consistent with a fiuid approximation up to second 
order. Our analysis is based on the careful definition of a locally Minkowski frame 
which thanks to the equivalence principle simplifies the treatment of local interactions. 
Though this method was already implicitly followed in the existing literature, we 
have performed a deep analysis of its implications for the gauge transformation, the 
SVT decomposition and the evolution equations. At first order it provides a more 
satisfactory understanding of the formalism used in the literature, whereas at second 
order it appears that it is completely necessary to master the formal aspects that 
we have presented in this paper in order to understand the physical meaning of the 
dynamical equations. For this purpose we have introduced a font and color based 
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notation which should clarify the formalism. The results that we have presented should 
now be the starting point to a numerical integration of the second order radiative 
transfer, but it appears already to be a huge and long term task if we want it to be 
computationally efficient. Similarly to the first order, approximate schemes should be 
first developed in order to obtain the main features, and this first step has already 
been taken in [44] for small angular scales. 
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Appendix A. Transformation rules for multipoles at second order in v 

Following the method of section 11.71 to compute the transformation rules of the 
radiation multipoles under a change of frame, we find the following rules up to second 
order in v for the multipoles used in section 16.51 (note that contrary to the choice we 
made in section [1.7.21 we express the result in function of the li"\p°) and not the 
/i"^ {p°) since this is what was required to transform the collision tensor in section [675)) 

W) = hip") - ^0(P°)n.v + + i/aP°)v.v + (n.v)^ 

+ Iv^hip") + Iv^I'dp") - i/rb°)«V.n - nip^v'^^.n 

+ Ih^ipny'v^ + In.ip^yv^ + , (A.l) 



%(/) = hjip") + -4{p°)v^,v,) - I(,{p'')v,) + I[.,{p°)v,) - I'(,{p'')v,)n.v 

- /^(p°)n.v - ^/,,(p°)v.v + ^4(p°)v.v + 3^/.';(p°)v.v + \l':^{p°) (n.v)^ 

AO 9n 9 
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+ ^E.^dPlv.^v'^ ^E',^^{p'')v,^v' ^E'^^,{p'')v,^v' 



k 



20 5 

+ l^E,,k{p")vk + i^E[^^{j>")vk 



63 ■'"^ ' 63 '^""^ ' 126 

+ lv,e'^\^B,)i{p'^) + \vke''\.^B'^)^{p^) (A.3) 

- {iBj)k,n{p )v VI - -e (^iBj^^^ip )v VI - — e ^iBj^f^„^{p )v vi . 

As for their energy integrated counterparts, they transform as 

4 2 2 ■ ■ 

2^0 = + -X0V.V - -ly + —lijv'v^ , (A.4) 

iij = X^j + lOX0U(iWj) - ^I{iVj) + Ik{iVj)v'^ , (A. 5) 

4- = - 3£k{,Vj^v'' + 2f.yv.v + 2«fce'='^,Sj), . (A.6) 

Appendix B. Sources terms in second order transformations 

The perturbation variables in the decomposition ()2.ip are extracted as foUows 

= - ^& (B.l) 

r(») - J_p ^o(") 

r(") _ p 



rr(n) _ J_p L ( <:! f:.J _ „(») 

-"mat — ^^^"M^fW I "a'Ol ^OklO I ff/j 



where n = 1, 2 is the order and where we have used the definitions 

P/ = A-la^ P/'' = <5^-' - A-l9^a•^ (B.2) 

Using this method we can read the source terms defined in equation (j2.3ip . which 
are quadratic in the gauge change variables T, L and the perturbation variables 

S^^T (T" + bUT' + {7i' + 2n^)T + + 2$') + diL'd^ (T - 2B - L') (B.3) 
+ T' (2r' + 4$) + diLd' (T' + HT + 2$) 

5** = - r (WT' + {n' + 27i;2)T - 2*' - m-^) - di {ht - 2*) d'l 

- i {S'-' - A-'d'd') Xij (B.4) 
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Sb = VjXi (B.5) 

Sb, = ^viXj (B.6) 

Se = (AA)-i - 1a5"^ A-.j (B.7) 



Se,, = 2P.^^P." (^,51.(5^ - ^5KL5'-'j Xij (B.8) 
= P.^P.M ('^I-'Jl' - ^'^A-L.J''') -^-/j , (B.9) 

with 

Xi = {T'di{2B + L' -T) + 2HTdi{2B + L' -T) 

+ 9 [2didjL + 2 (HT - 2^-) Su + 4 {Hu + 9/9,/ S)] 
+ d 'diLdj {2B + L' -T) + d'' Ldjdi {2B + L' - T) 

+a/r(-4$ - 2T' - 2nT) + Tdi{2B' + L" - T')} , (B.IO) 

Xu = {dj {2B + L' -T) diT + Tdidj{L' + 2HL) 

+ did^L [2dKdjL + AdKdjE + AHkj + {2'HT - A^)5kj] 
+ T {2H'jj + 2didjE' + AUHu + AUdidjE) 

+d'^LdK {didjL + 2Hij + 2did.,E)] . (B.ll) 

As for the matter perturbation variables, the source terms in the transformation rules 
are 

Sp ^T{p"T + p'T' + 25p')+d'Ldi{25p + p'T) (B.12) 

Sp ^T{P"T + P'T' + 26P')+d'Ldi{26P + P'T) (B.13) 
Sv = Psi [HTd'{L' - 2V) + Td\2V' - L") + d \L' - 2V)djd' L 

+L-'djd'{2V - L') + d'L' {HT + T' + 2$)] (B.14) 
Sy„ = P^f [HTd'iL' - 2V) + Td\2V' - L") + d\L' - 2V)djd'L 

+L-'djd\2V - L') + d'L' {HT + T' + 2$)] (B.15) 

S^ij = 2T (tt^')' + 2d'^'LdKn'' - 2n"'^dKd-'L - 2t:-'^ OkO'L . (B.16) 

Appendix C. Perturbation of tetrads and Ricci rotation coefficients 

Appendix C.l. The perturbation of tetrads 

At first order the coefficients of Rab and Sat are (remembering that we discard the 
first order vector modes) 

xR^;j = xSi'J = $(1) (C.f) 
xR^'> = - xS^ ^ 

xR^=-S^=¥^^5,,-d,d,Ei^^H^ 
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Wc can read directly from these expressions the transformation rules for the tetrad 
when going from a gauge X to a gauge Y 



ye, - J 



,(1) 



(1)1 



At second order the coefficients of Rat and Sab are 
Ril^ = - 3$2 + diBd^B 

- - djB^^'> - Bf^ + (2$ - m)diB + Ad'B {djdjE + Hi,,) 



X ^ ^oo 
,(2) 



(C.2) 



(C.3) 



X 



X 



.(2) 
Ho 

4^^ 



X 



(2) 

io 

.(2) 
'ifc 







j(2) 



X ^oo 

■xSi 



+ 3 + ijf ) (^L^A-^ + i/LA-) - 6* {didKE + Hik) 

$(2) _ $2 ^ diBd^B 

- diB^^'> - 2^diB + 2a'^B (a/^ji; + Hij) 



The transformations rules for the tetrads can then be read 



(.ve?') 



(.vef) 



(C.4) 



'-a^r(B(2)) _ r(B^(2)) ^ [2r($) - 4r(i')] a^T(s) 
+4a''r(B) [a'ajr(£;) + h[j] } e, 

+3 [97a''r(i;) + nj] [d"djT{E) + h'}] } e^. (c.s) 

Appendix C.2. The perturbation of Ricci rotation coefficients 

If we use a tetrad basis, then the covariant derivative is characterized by the Ricci 
rotation coefficients [75] defined by 



^abc 



l^Me^e^'V^e^ 



(C.6) 



They are related to the Christoffel symbols according to 



pC A b 



dee' 



We choose a background tetrad field which is adapted to the Cartesian coordinate 
system (see section [3T3| . We then expand the Ricci rotation coefficients by expanding 
both the Christoffel symbols and the tetrads. Since the Ricci rotation coefficients are 
antisymmetric in their two last indices. lOooo = ^^ioo = up to any order. For the 
background metric the only non- vanishing components are 



= -LOi.o = ()• 



(C.S) 



We can check that at this order the vectors in the tetrad commute (since they arise 
from a coordinate system) as 



(C.9) 
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At first order, restricting to the Newtonian gauge and neglecting the vector 
perturbations, the non- vanishing components are 



1 



a 

(1) 2 



UJ 



lOJ 

(1) 



jik 



jki 



a 



^ = 0. 



(C.IO) 

As for the second order, the components that we have used in this paper are in the 
Newtonian gauge 



LO 





(2) 

— CJ ■ = 


1 - 


OOl 


OlO 


a - 













7i;$(2) + ,5(2)' 



(C.ll) 



Si J - 3H<^^Sij 



+2[H'jj, - ^'5jk] [Hf - mf] } . (C.12) 
It can be checked that the vectors in the tetrad do not commute at the perturbed level 



as ^ ^c6a and w),;^^ ^ 



,(1) 



,(2) 



Appendix D. Useful formulas for extracting the normal modes 

From the general composition rule (|7.29p . we obtain the following useful particular 
cases which can be used to extract the normal modes on quadratic terms 

-'^tyt^yn.h — r— ; iViAiJ^ft , 



^+1 



A .,*M^;(^+i)(m±i) _ v/(£+l±m)(£ + 2±?7t) 



(2^-1) 



b 1 



iv^-i V2{2e-i) ' 



(D.l) 
(D.2) 
(D.3) 
(D.4) 



In order to extract normal modes in quadratic terms involving a curl, we also need 
the following useful formulas 



; '\jtm '\jllra 



^(±1) 6c y(mTi) ^ , 1 / (^ + mTl + l)(^-(mTl)) ,,^(^) 

'6 (ai -'^ a(^_i) )e 2 ^ ' 



(D.5) 
(D.6) 
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